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ABSTRACT 


A  3— year  (1986-89)  micromechanical  research  at  Rensselaer  Polytechnic 
Institute  (RPI)  on  the  behavior  and  modelling  of  granular  media  is  summarized. 
The  final  objective  is  to  develop  a  constitutive  law  for  granular  soil  based  on  the 
particulate  nature  of  the  material. 

This  is  accomplished  by  a  systematic,  mostly  analytical  approach  to  the 
problem,  starting  from  the  response  of  the  contact  between  two  elastic  rough 
spheres  subjected  to  arbitrary  normal  and  tangential  forces,  and  continuing  with  the 
response  of  regular  and  random  arrays  of  spheres.  The  following  tasks  were 
completed:  a)  study  and  compilation  of  the  differential  stress— strain  relationships  of 
several  regular  arrays  of  identical  quartz  spheres;  b)  use  of  the  Self  Consistent  and 
Nonlinear  Finite  Element  methods  to  calculate  the  small  strain,  monotonic  and 
cyclic  stress— strain  behavior  of  random/regular  arrays  of  identical  quartz  spheres 
loaded  isotropically  and  anisotropically,  including  wave  velocity  predictions;  and  c) 
use  of  Nonlinear  Distinct  Element  simulations  of  two-dimensional  random  arrays  of 
quartz  spheres  to  determine  their  small  and  large  strain  response,  including  the 
initial  position  and  subsequent  translation  and  distortion  of  yield  surfaces. 

Throughout  the  research,  comparisons  were  made  between  analytical 
predictions  and  experimental  measurements  on  sands  and  other  granular  media 
reported  in  the  literature.  In  addition,  several  computer  controlled,  hollow- 
cylinder,  monotonic  and  cyclic  axial— torsional  tests  were  conducted  on  glass  bead 
specimens,  along  stress  paths  similar  to  those  used  in  the  Distinct  Element 
simulations,  to  verify  and  supplement  the  results  of  the  calculations. 
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Finally,  a  constitutive  iaw  for  granular  media  is  proposed  and  its  main 
features  are  discussed.  The  law  is  basically  the  stress— strain  equivalent  of  the 
force— deformation  model  for  the  contact  between  two  spheres  developed  at  RPI, 
enhanced  to  incorporate  dilation  and  the  distortion  of  yield  surfaces  due  to  prestrain 
which  has  been  observed  in  granular  media.  The  proposed  stress— strain  law 
includes  an  infinite  number  of  yield  surfaces  of  conical  shape  initial  parallel  to  each 
other,  as  well  as  modified  normality  and  kinematic  strain  hardening  rules. 
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and  (d)  Textured  Magnesium  (Kelley  and  Hosford,  1968). 

Figure  60.  Initial  and  Subsequent  Yield  Surfaces  of  the  531— Sphere  Medium  in 
the  TVh  1/2  (<Tv-<Th)  Space  calculated  Using  Distinct  Element 
Program  CONBAL-2.  Yield  Surfaces  Defined  as  the  Locus  of  all 
Points  with  70Ct  =  0.05%. 

Figure  61.  Conical  Y  ield  Surfaces  and  Current  r— plane  for  Constitutive  Law 
of  Granular  Soil  of  Eqs.  1—7.  Position  of  Yield  Cones  Correspond 
to  Loading  from  O  to  A  followed  by  Stress  increment  AB.  Point  B 
is  on  the  T=plane.  The  Value  of  p  at  the  Apex  of  Any  Cone  (e.g., 
A,  B)  is  pi. 
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Figure  65. 


Schematic  Representation  of  (a)  Translation  and  (b)  Distortion  of 
a  Yield  Surface  (Yen  1979). 

Predicted  and  Measured  Yield  Surfaces  in  Aluminum  (Yen  and 
Eisenberg  1987). 

Yield  Surface  in  principal  Stress  Space.  Note  the  Rotation  of  the 
Surface  Around  its  Apex  (Prevost  1985). 

Proposed  Modification  of  the  Hardening  Rule  to  Allow  Limited 
Rotation  Around  the  Yield  Surface  Apex  in  Order  to  Account  for 
Dilation. 
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1.  INTRODUCTION 


This  final  report  presents  a  summary  of  the  work  conducted  by  the  authors  in 
a  three  year  (1986-89)  AFOSR  sponsored  micromechanical  research  at  RPI  on  the 
behavior  and  modelling  of  granular  media.  Four  progress  reports  were  submitted  to 
AFOSR  for  years  1986—87  and  1987—88  and  are  included  in  the  list  of  references 
(Petrakis  and  Dobry  1987,  1987a,  1988,  and  Petrakis,  Dobry  and  Ng  1988).  The 
final  objective  of  the  research  is  to  develop  a  constitutive  model  for  granular  soils 
expressed  in  traditional  continuum  mechanics  terms,  but  grounded  on  the 
particulate  micromechanical  behavior  of  the  material.  As  sketched  in  the  flow  chart 
of  Table  1,  this  is  accomplished  by  a  systematic,  mostly  analytical,  approach  to  the 
problem,  starting  from  the  response  of  the  contact  between  two  elastic,  rough 
spheres  subjected  to  arbitrary  normal  and  tangential  forces,  and  continuing  with  the 
response  of  regular  and  random  arrays  of  spheres.  Hollow-cylindrical  axial- 
torsional  tests  on  glass  beads  are  also  performed  to  supplement  the  analytical 
calculations.  These  analytical  and  experimental  results  are  then  used  to  formulate 
the  desired  constitutive  model  for  granular  media. 

The  problem  of  the  contact  between  two  identical  spheres,  originally  solved  by 
Mindlin  and  Deresiewicz  for  some  special  cases,  had  been  solved  at  RPI  in  1984  for 
arbitrary  loading  through  an  incremental  elastic-plastic  model.  This  general 
solution  was  then  approximately  adapted  for  unequal  spheres,  and  implemented  and 
used  in  this  research  to  study  a  number  of  regular  and  random  arrays  of  quartz 
spheres  utilizing  several  analytical  and  numerical  techniques.  These  techniques 
include  the  pressure  dependent  Self-Consistent  method,  the  nonlinear  Finite 
Element  Method,  and  the  nonlinear  Distinct  Element  Method.  The  results  obtained 
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ranged  from  small  to  large  strains  and  included  both  monotonic  and  cyclic  loading 
response. 

Definite  results  and  conclusions  are  presented  in  this  report  for  most  of  the 
tasks  included  in  Table  1.  This  includes  the  two— sphere  problem,  the  study  of 
regular  arrays,  and  the  use  of  the  Self— Consistent  and  Finite  Element  methods  to 
study  random  arrays  of  spheres. 

The  purpose  of  the  experimental  effort  was  to  verify  the  analytical  and 
numerical  predictions.  This  experimental  task,  also  included  in  Table  1,  consisted 
of  eight  hollow— cylinder  axial  torsional  tests  on  glass  beads  to  observe  the  shape  and 
evolution  of  yield  surfaces  of  granular  media  in  stress  space.  Distinct  Element 
numerical  simulations  on  random  arrays  of  spheres  were  also  conducted  along  the 
same  stress  paths  used  in  the  experiments.  Both  experiments  and  simulations 
suggest  that  the  yield  surface(s)  of  a  granular  medium  distort(s)  in  the  direction  of 
loading  but  not  in  other  directions,  and  this  finding  is  incorporated  into  the 
constitutive  model  proposed  at  the  end  of  the  report. 


1.1  Brief  Review  of  Constitutive  Relations  in  Soils 

Over  the  past  30  years  considerable  attention  has  been  given  to  the 
development  of  constitutive  laws  for  engineering  materials  (Hill  1950,  Prager  1955, 
Mroz  1967,  Dafalias  and  Popov  1976,  Drucker  and  Palgen  1981).  Among  other 
formulations  the  existing  models  are  based  on  the  theories  of  elasticity, 
hypoelasticity,  plasticity  and  viscoplasticity.  Despite  the  large  number  of  models, 
there  is  no  consensus  yet  within  the  research  community  on  the  best  approach. 
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However,  the  models  based  on  the  theory  of  plasticity  or  viscoplasticity  appear  to  be 
most  promising. 

Most  of  the  proposed  models  for  soils  are  based  on  the  incremental  theory  of 
plasticity.  In  these,  the  total  strain  increment  is  equal  to  the  sum  of  elastic  and 
plastic  strain  increments,  dry  =  defj  +  dr?j,  with  these  increments  being  rate 
independent  (Drucker  and  Prager  1952,  Reyes  1966,  Chen  1975,  Lade  and  Duncan 
1975,  Prevost  1978,  Hardin  1978). 

A  variety  of  associative  and  non-associative  flow  rules  have  been  proposed  for 
the  plastic  strain  increment,  of  the  form: 

<*'*i  -  dA  M 

where  dA  is  a  coefficient  of  proportionality  and  g(cnj)  is  the  plastic  potential 
function,  which  may  or  may  not  coincide  with  the  yield  function,  f(cxij)  at  which 
plastic  strains  develop. 

In  the  simplest  type  of  elastic— plastic  model,  there  is  only  one  yield  surface. 
States  of  stress  below  that  surface  are  assumed  to  be  elastic.  However,  soils  develop 
plastic  strains  even  at  very  small  strains;  to  allow  for  this  behavior,  a  wide  variety 
of  strain-hardening  laws  have  been  proposed,  including  families  of  yield  surfaces 
and  specific  strain  hardening  yield  rules.  In  some  of  these  models  the  elastic  region 
is  completely  eliminated,  allowing  for  plastic  flow  at  very  low  strain  levels  (Prevost 
1978). 

An  important  aspect  of  the  development  of  elastic— plastic  models  is  the 
definition  of  the  strain  hardening  law,  which  defines  the  modifications  in  the  yield 
surface(s)  in  the  course  of  plastic  flow.  This  is  critical  for  reversible  loading,  where 


4 


the  type  of  strain  hardening  determines  the  stress-strain  behavior  after  load 
reversals.  In  some  of  the  older  general  plasticity  models  developed  for  monotonic 
loading,  isotropic  hardening  is  assumed  (Hill  1950),  with  the  yield  surface(s) 
expanding  in  size  as  the  stresses  increase.  With  isotropic  hardening,  a  large  amount 
of  load  reversal  is  required  for  additional  yielding  to  occur,  in  contradiction  with  the 
observed  behavior  in  the  laboratory. 

A  better  alternative  is  provided  by  the  kinematic  hardening  law  proposed  by 
Ishlinsky  (1954)  and  Prager  (1955),  which  assumes  that  the  yield  surface  translates 
in  stress  space  without  changing  shape  or  size  during  flow.  Models  with  kinematic 
strain  hardening  are  in  better  agreement  with  experimental  results  than  models  with 
isotropic  hardening. 

Some  of  the  earliest  and  most  popular  formulations  for  sand  have  been  based 
on  the  "cap  model"  of  DiMaggio  and  Sandler  (1971).  As  versatile  as  "cap  models" 
may  be,  they  have  not  been  successful  in  accurately  modelling  cyclic  loading 
conditions. 

Similar  limitations  apply  to  other  models  (Bala.  •  Rohani  1979),  and  can 
be  argued  that  the  existing  plasticity  models  of  the  "cap"  type  are  only  adequate  for 
monotonic  loading  of  isotropic  soil.  In  an  effort  to  overcome  this,  a  variety  of 
constitutive  laws  have  been  proposed  which  incorporate  a  combination  of  isotropic 
and  kinematic  hardening  (Mroz  1967,  Lade  1977).  Although  these  more  recent 
theories  represent  a  considerable  advancement  over  the  "cap"  models,  they  too  have 
drawbacks.  These  include  the  use  of  "a  priori"  geometrical  hardening  rules,  and  the 
fact  that  they  do  not  predict  either  the  inherent,  elastic  anisotropy  of  the  soil  or  the 
stiffening  effect  measured  after  many  load  reversals. 
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This  elastic  (inherent)  anisotropy,  which  is  most  significant  for  anisotropically 
consolidated  sand,  has  been  measured  in  sand  by  Stokoe  and  his  coworkers  (Knox  et 
al.  1982,  Koppermann  et  al.  1982),  while  Dafalias  (1979)  has  discussed  its 
implications  for  modelling.  In  their  experiments  Stokoe  and  his  coworkers  proved 
that  the  wave  propagation  velocity  depends  on  the  principal  stresses  in  the  direction 
of  wave  propagation  and  particle  motion  (Fig.  1)  and  not  on  the  mean  stress  as 
proposed  by  Hardin  and  Black  (1964)  and  Seed  and  Idriss  (1970).  This  elastic 
anisotropy  has  not,  to  the  best  of  the  authors’  knowledge,  yet  been  incorporated  into 
any  plasticity  model.  The  existing  plasticity  models  also  do  not  predict  the 
dramatic  stiffening  occurring  in  a  granular  medium  which  has  been  precycled  for 
thousands  or  millions  of  cycles  at  small  strains  (7  ~  6  X  10"4  )  without  experiencing 
significant  volumetric  changes  (Fig.  2,  see  also  Drnevich  1967).  These  two  examples 
are  the  result  of  micromechanical  phenomena  and  are  best  modelled  once  these 
micromechanical  phenomena  are  understood. 

Another  way  to  say  the  same  thing  is  to  state  that  until  now  plasticity  models 
for  soil  have  been  mostly  phenomenological.  They  have  been  typically  developed 
from  a  manageable  mathematical  formulation,  and  they  have  been  calibrated  and 
modified  by  interpreting  macroscopic  experimental  results.  The  underlying 
micromechanical  phenomena  have  not  been  systematically  considered.  As  a  result, 
the  existing  plasticity  models  for  soils  are  in  need  of  constant  refinement  when 
needed  for  cases  very  different  from  the  one  the  model  was  originally  developed  and 
calibrated  for. 

The  current  situation  in  metal  plasticity  is  quite  different.  Although 
modelling  of  the  nonlinear  behavior  of  metals  started  on  a  similar  phenomenological 
basis,  there  has  been  a  shift  in  the  last  20  years  or  so  toward  formulating  the  metal 
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response  with  due  consideration  of  micromechanical  principles  (Budiansky  and  Wu 
1962,  Lin  and  Ito  1965,  1966).  Recently  this  has  been  enhanced  by  specific 
experiments  and  micromechanical  (electron  microscopy)  measurements  (Stout  et  al. 
1985,  Helling  et  al.  1986).  The  situation  is  analogous  in  the  modelling  of  more 
complex  composite  materials,  where  experiments  and  micromechanical  analytical 
simulations  are  combined  to  create  the  corresponding  constitutive  law  (Dvorak 
1987,  Dvorak  et  al.  1988). 

Although  metal  properties  are  not  pressure  dependent  and  their  material 
symmetry  does  not  change  as  much  with  loading  as  in  soils,  the  stress— strain 
behavior  of  metals  and  soils  is  similar  in  several  respects.  As  a  result,  most  current 
soil  plasticity  models  are  modified  versions  of  popular  phenomenological  metal 
models.  Notable  examples  include  the  Mroz  (1967)  model  for  metals  and  the 
Prevost  (1978)  model  for  undrained  loading  of  clay,  as  well  as  the  bounding  surface 
model  used  by  Dafalias  and  Herrmann  (1982).  Unfortunately,  no  plasticity  model 
exists  for  soil  resulting  from  the  combination  of  specific  laboratory  experiments  and 
micromechanical  principles  including  numerical  simulations  of  the  behavior  of 
granular  arrays  under  load. 

Recently,  the  geotechnical  groups  at  the  University  of  Colorado  (Klisinski  et 
al.  1988),  US  Army  WES  (Peters  1988),  and  University  College  of  London  (Arthur 
et  al.  1988),  working  together  in  a  concentrated  effort,  developed  a  constitutive  law 
which  is  based  on  a  series  of  innovative  3-D  laboratory  experiments  performed  for 
this  purpose  (Alawi  et  al.  1988).  This  may  be  the  first  time  that  a  comprehensive 
experimental  investigation  is  performed  in  Soil  Mechanics  which  attempts  to  verify 
widely  accepted  rules  and  assumptions  related  to  hardening,,  yield  surface  shape 
during  loading,  normality  rule,  etc  (see  Fig.  3).  The  resulting  model  (Klisinski  et  al. 
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1988)  is  based  on  the  bounding  surface  model  of  Dafalias  and  Popov  (1976)  and  on 
the  theory  of  "fuzzy"  sets.  The  model  was  used  with  moderate  success  to  predict 
the  behavior  of  sand  in  the  International  Workshop  on  Constitutive  Laws  at  Case 
Western  Reserve  University  (Saada  1987). 

In  the  above  approach  many  lab  specimens  were  used  to  define  one  yield 
surface,  which  is  satisfactory  for  monotonic  loading  but  omits  the  effect  of 
prestraining  which  in  metals  is  known  to  play  an  important  role  in  cyclic  loading. 
Moreover,  the  development  of  the  corresponding  law  is  again  phenomenological  and 
includes  no  micromechanical  considerations,  relying  instead  on  the  interpretation  of 
the  macroscopic  results  of  the  above  innovative  experiments. 


1.2  Micromechanical  Interpretation  of  the  Deformation  Mechanism  in 
Polvcrvstalline  Aggregates. 

For  some  purposes,  granular  soil  can  be  modelled  as  a  pressure— dependent 
polycrystalline  aggregate.  Regular  arrays  of  identical  spheres  studied  by 
Deresiewicz  (1958,  1958a)  and  Petrakis  and  Dobry  (1989)  behave  like  pressure- 
dependent  single  crystals.  Random  arrangements  of  these  regular  arrays 
(poly crystals)  have  been  successfully  used  to  simulate  the  small  strain  behavior  of 
sand  by  means  of  the  Self  Consistent  and  Finite  Element  techniques  (Petrakis  and 
Dobry  1986,  1987).  Thus  it  is  useful,  from  the  viewpoint  of  the  development  of  a 
constitutive  law  for  granular  soil,  to  review  some  findings  on  poly  crystalline 
aggregates  such  as  metals. 
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Starting  in  the  1950’s,  researchers  have  simulated  the  elastic— plastic, 
stress-strain  behavior  of  such  polycrystalline  aggregates  through  analytical, 
semi— analytical,  and  numerical  micro  mechanical  techniques  (Hershey  1954, 
Budiansky  and  Wu  1962,  Lin  and  Ito  1965,  1966,  Hill  1967,  Canova  et  al.  1985). 
This  micromechanical  work  has  not  supported  the  continuum  mechanics  hypotheses 
of  either  pure  kinematic  or  isotropic  strain  hardening  hardening  behavior,  but  has 
predicted  instead  a  combined  translation  (kinematic  hardening)  and  distortion  of 
the  yield  surfaces  in  the  direction  of  loading.  This  micromechanical  prediction  has 
been  verified  by  several  experiments  (Naghdi  et  al.  1958,  Phillips  1968,  Phillips  and 
Tang  1972,  Phillips  et  al.  1974). 

The  micromechanical  approach  commonly  used  to  analyze  the  elastic— plastic 
behavior  of  polycrystals  assumes  that  they  are  an  assemblage  of  equal  anisotropic 
monocrystals  (single  crystals),  randomly  oriented  in  space  (Fig.  4a).  This  results  in 
an  isotropic  polycrystal  if  the  spatial  distribution  of  the  orientations  is  statistically 
uniform.  A  monocrystal  has  n  sliding  planesv  with  each  plane  having  m  sliding 
directions,  and  with  every  sliding  direction  corresponding  to  a  pair  of  parallel  yield 
planes  in  stress  space.  In  the  limiting  case  in  which  an  infinite  number  of  possible 
crystal  orientations  is  assumed,  this  infinitely  sided  polyhedron  becomes  a  curved 
yield  surface. 


(*)  In  the  literature  dealing  with  polycrstals,  they  are  called  "slip  planes."  Here 
they  will  be  referred  to  as  "sliding  planes,"  as  we  will  reserve  the  use  of  the 
word  "slip"  to  denote  a  different  phenomenon  when  discussing  the  contact 
between  two  elastic  rough  spheres 
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Plastic  strain  in  the  aggregate  is  caused  by  sliding  of  one  of  the  sliding  planes 
occuring  in  a  family  of  similarly  oriented  crystals.  After  sliding  has  occurred  in  a 
number  of  these  families,  each  surface  of  the  polyhedron  mentioned  above  expands 
and  shifts  differently.  These  sliding  directions  are  all  more  or  less  parallel  to  the 
direction  of  the  plane  of  the  maximum  shear  stress  acting  on  the  aggregate.  As  the 
aggregate  is  loaded  further  beyond  the  elastic  range,  more  crystals  and  crystal 
families  slide,  and  increasingly  more  yield  planes  pass  through  the  loading  point  on 
the  yield  surface.  These  yield  planes  of  different  orientations  intersect  at  that  point 
on  the  yic’d  surface  and  form  a  corner  or  vertex  (Fig  4b). 

This  vertex,  which  is  particularly  important  for  stress-strain  modelling  during 
cyclic  loading,  is  not  easily  observed  during  testing.  One  reason  is  that  the  very 
large  number  of  monocrystal  orientations  smoothes  the  effect,  which  appears  as  a 
"smooth  vertex"  or  distortion  of  the  yield  surface  in  the  direction  of  loading,  rather 
than  a  sharp  corner.  Most  important,  this  "vertex"  does  not  appear  at  all  if  purely 
monotonic  tests  are  performed,  but  it  shows  up  once  the  loading  is  reversed  at  least 
once.  This  distortion  of  the  yield  surface  associated  with  the  vertex  reflects  the 
"memory"  the  material  has  of  prestraining  in  the  direction  of  loading.  The 
existence  of  this  yield  surface  distortion  in  metals  has  been  observed  experimentally 
by  a  number  of  researchers  in  several  polycrystalline  aggregates,  including 
aluminum,  aluminum  alloys,  brass  and  magnesium  (Naghdi  et  al.  1958,  Phillips 
1968,  Kelley  and  Hosford  1968,  Phillips  et  al.  1970,  Shiratori  et  al.  1976,  Helling 
et  al.  1986). 

The  late  Professor  Phillips  developed  a  testing  procedure  to  seek  and  compute 
the  initial  and  subsequent  yield  surfaces  of  aluminum  and  their  motion  in  stress 
space.  This  procedure  is  now  widely  used  in  experimental  plasticity  studies  of 
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metals  and  metal  matrix  composites  (Rousset  1985,  Stout  et  al.  1985,  Dvorak  1987, 
Dvorak  et  al.  1988).  The  experiments  are  typically  conducted  by  applying  a 
combination  of  tensile  (a)  and  torsional  shear  (r)  stresses  to  a  hollow  cylindrical 
specimen,  in  a  sequence  similar  to  that  shown  in  Fig.  5.  In  these  tests,  the  loading 
stops  and  reverses  as  soon  as  a  point  on  the  yield  surface  is  reached  —  so  that  the 
entire  yield  surface  may  be  determined  —  as  defined  by  a  certain  deviation  from  the 
linear  portion  of  the  stress-strain  curve.  Fig.  6  clearly  shows  the  characteristic 
distortion  of  the  initial  yield  surface  in  the  direction  of  loading.  While  the  yield 
surface  (for  a  given  temperature)  in  the  r—  a  space  is  an  ellipse,  the  subsequent  yield 
surfaces  have  distorted  and  become  pointed  in  the  direction  of  loading  (a),  while 
becoming  flatter  in  the  opposite  direction  (b).  As  a  result,  the  size  of  the  yield 
surface  shrinks  in  the  direction  of  loading  while  staying  constant  in  the  other 
direction. 

Laboratory  results  such  as  these  have  made  possible  the  linking  of  the 
micromechanical  theory  with  experiments  (Stout  et  al.  1985,  Helling  et  al.  1986), 
and  have  led  to  a  new  family  of  constitutive  laws  (Phillips  and  Weng  1975, 
Eisenberg  and  Yen  1981,  1984,  Yen  and  Eisenberg  1987)  which  incorporate  the 
above  findings. 


1.3  Micrcmechanical  Behavior  of  Granular  Soil 

As  mentioned  earlier,  the  behavior  of  a  sand  aggregate  is  similar  to  that  of  a 
polycrystal,  since  the  individual  groups  of  grains  or  grain  packings  within  the  sand 
may  be  considered  in  first  approximation  to  behave  like  randomly  oriented  crystals. 
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The  main  difference  is  that  the  properties  of  these  packings  are  now  pressure 
dependent.  For  example,  a  simple  cubic  array  of  equal  spheres  (Fig.  7)  is  a 
pressure— dependent  monocrystal  having  three  sliding  planes  (n=3),  with  each 
sliding  plane  containing  two  sliding  directions  90°  apart  (m=2). 

As  in  the  case  of  polycrystalline  aggregates,  each  sliding  plane  in  each  of  the 
packings  corresponds  to  a  pair  of  parallel  yield  planes  in  stress  space.  The 
macroscopic  yield  surface  of  the  material  is  the  surface  bounding  the  yield 
polyhedron,  the  sides  of  which  are  formed  by  the  intersection  of  these  yield  planes. 
Plastic  strain,  defined  as  irrecoverable  deformation,  is  the  result  of  a  slide  in  at  least 
one  of  these  packings.  A  conjecture  to  be  discussed  later  in  this  report  is  that  the 
yield  surface  of  a  granular  medium  distorts  during  loading,  and  forms  a  vertex  in 
the  direction  of  loading  through  a  mechanism  similar  to  those  in  metals. 

Granular  media,  despite  this  similarity  with  metals,  are  also  significantly 
different  due  to  the  pressure-dependent  properties  of  the  packing.  Specifically,  the 
amount  of  slide  in  each  packing  depends  on  the  normal  pressure  acting  on  the  plane, 
as  do  the  moduli  and  the  symmetry  of  the  material.  Moreover,  granular  media 
experience  dilation  under  shear  which  is  not  necessarily  present  in  polycrystalline 
aggregates.  Finally,  soils  exhibit  nonlinear  inelastic  stress-strain  behavior  even  at 
small  strains.  Therefore,  strictly  speaking,  a  cohesionless  aggregate  does  not  have  a 
clear  linear  elastic  region  followed  by  a  region  in  which  permanent  deformations 
occur,  like  the  aluminum  tested  by  Phillips.  However,  granular  soils  experience 
little  or  not  sliding  (gross  sliding)  between  particle  contacts,  and  thus  exhibit 
nondestructive  behavior  and  no  dilation  up  to  the  so  called  "threshold  shear  strain," 
7t  ~  10 (Dobry  et  al.  1982).  Therefore,  in  granular  soils,  loading  below  7t  has 
some  important  features  in  common  with  elastic  loading  in  metals  before  the  initial 
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yield  surface  is  reached.  However,  in  soils  the  loading  below  7t,  though 
nondestructive,  is  inelastic  and  does  include  some  plastic  yielding  due  to  localized 
slipping  within  the  intergranular  contact  areas.  If  this  localized  slipping  effect 
below  7t  is  neglected  in  first  approximation,  one  possible  definition  of  an  initial 
yield  surface  in  a  cohesionless  soils  could  be  the  surface  in  stress  space  where  7t  is 
reached,  there  is  significant  sliding  between  particles,  and  permanent  volumetric 
deformation  starts  occurring  due  to  shear.  This  may  be  measured  by  monitoring  the 
volume  in  drained  experiments  or  the  residual  pore  pressure  in  undrained 
experiments.  These  inelastic  deformations  are  typically  the  result  of  dilatancy  in 
the  granular  medium  and  would  not  occur  without  this  sliding  between  particles. 
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2.  THE  MICROMECHANICAL  APPROACH 

2.1  Problem  of  Contact  Between  Two  Elastic  Rough  Spheres 

The  problem  of  the  contact  between  two  elastic  rough  spheres  is  of  great 
importance  to  this  research.  The  small  strain  behavior  of  a  granular  medium 
formed  by  spherical  particles  is  controlled  jointly  by  the  elasticity  and  the  friction 
coefficient  of  the  contacts,  with  friction  becoming  more  important  at  very  large 
shear  strains  approaching  failure  of  the  medium.  Therefore,  a  brief  review  of  the 
contact  problem  will  be  presented  herein,  as  well  as  some  aspects  of  the  recent 
numerical  solution  for  the  contact  problem  developed  at  RPI. 

The  problem  of  contact  of  two  elastic,  elliptical,  semi-infinite  bodies  subjected 
to  a  normal  force  was  first  studied  by  Hertz  (1882),  with  this  solution  including  as  a 
special  case  the  normal  compression  of  two  spheres.  Hertz  demonstrated  for  the 
first  time  that  the  normal  force— deformation  behavior  at  the  contact  is  nonlinear 
elastic.  Subsequently,  all  work  on  the  same  topic  was  concerned  with  the  loading  of 
bodies  by  normal  forces,  until  Cattaneo  (1938),  Mindlin  (1949),  and  Mindlin  and 
Deresiewicz  (1953)  addressed  the  problem  of  the  contact  of  two  identical  elastic, 
rough  spheres  subjected  to  a  combination  of  normal  and  tangential  forces,  and 
presented  a  number  of  closed  form  solutions  for  each  of  several  specific  load 
histories.  Walton  (1978)  studied  the  problem  of  the  oblique  compression  of  two 
elastic,  rough  spheres  when  both  tangential  and  normal  forces  are  applied 
simultaneously.  More  recently,  Szalwinski  (1985)  addresed  the  problem  of  the 
contact  of  two  identical  elastic  bodies  which  form  an  elliptical  contact  area,  and 
compared  his  solution  to  those  in  the  literature. 
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The  general  case  of  two  identical,  elastic,  rough  spheres  subjected  to  a  normal 
force  N  followed  by  a  tangential  force  T  (Fig.  8),  solved  by  Cattaneo  (1938)  and 
Mindlin  (1949),  is  a  problem  of  the  linear  theory  of  elasticity.  Since  the  solution 
yields  an  infinite  shear  stress  at  the  edge  of  the  contact  area,  a  slip  needs  to  be 
prescribed  at  this  edge,  which  transforms  the  formulation  into  a  mixed  boundary 
value  problem  where  contact  stresses  and  displacements  are  prescribed.  This 
permanent  set  produced  by  the  slip  induces  a  nonlinear  behavior,  which  is  different 
from  that  computed  by  Hertz,  since  it  is  accompanied  by  energy  dissipation.  As 
demonstrated  by  Mindlin  and  Deresiewicz  (1953),  due  to  this  slip  now  the 
force-deformation  relation  depends  on  the  entire  past  history  of  the  loading  as  well 
as  on  the  instantaneous  rates  of  change  of  the  normal  and  tangential  forces.  A 
typical  force-deformation  curve  for  two  spheres  in  contact  under  a  constant  normal 
force,  N,  subjected  to  a  monotonicaily  increasing  tangential  force,  T,  is  shown  in 
Fig.  8b,  where  the  nonlinear,  yielding  behavior  can  be  clearly  observed. 

All  of  the  above  suggest  that  a  phenomenological  plasticity  model  could 
describe  this  nonlinear  behavior.  Such  a  formulation  would  provide  the  long 
awaited  (Deresiewicz  1958a)  "general'1  solution  to  the  problem  of  the  contact  of  two 
elastic,  rough  spheres  subjected  to  an  arbitrary  force  history,  which  in  turn  could  be 
used  in  numerical  simulations.  This  has  been  achieved  recently  at  RPI  through  a 
constitutive  law  (Seridi  and  Dobry  1984,  Dobry  et  al.  1989)  for  the  force- 
deformation  behavior  of  two  identical  elastic,  rough  spheres  in  contact  under  a 
combination  of  arbitrarily  varying  normal  and  tangential  forces,  which  was 
implemented  through  program  CONTACT  (Dobry  et  al.  1989).  This  model  is  based 
on  the  incremental  theory  of  plasticity,  uses  an  infinite  number  of  yield  surfaces  and 
assumes  kinematic  hardening.  Therefore,  its  main  features  are  very  similar  to  those 
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of  the  plasticity  stress— strain  models  for  engineering  materials  described  previously. 
These  features  of  the  RPI  contact  model  are  presented  in  detail  by  Seridi  and  Dobry 
(1984)  and  Dobry  et  al.  (1989),  and  are  summarized  below: 

Yield  Condition: 


(T*  —  x)2+  (Ty  —  y)!=  f2(N  —  Ni)a  (2) 

where  x  =  x  i  +  y  j  +  Nikis  the  position  vector  of  the  apex  of  the  conic  surface, 
and  P  =  Txi  +  Tyj  +  Nk  is  the  current  force  point  (Fig.  9).  Since  the  contact  of 
the  two  spheres  subjected  to  an  increasing  tangential  force  under  constant  N  is 
continuously  slipping,  there  must  be  an  infinite  number  of  yield  cones.  The  elastic 
region  is  in  the  inside  of  the  cone  w^ose  apex  is  at  the  current  force  point.  The 
failure  surface  is  the  outer  cone,  defined  as  follows: 

T2=  T|  +  Ty  =  f2N2  (3) 


Flow  Rule: 


For  a  given  force  increment  dP  =  d$  —  d'T,  the  increment  of  displacement 
between  the  centers  of  the  two  spheres  is : 


dO  =  ddxi  +  dtfyj  +  dak  =  d^  +  dale 


(4) 
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The  value  of  dC  is  given  by: 


where:  n=  normal  unit  vector  to  the  yield  circle  at  the  current  force  point. 

t=  tangential  unit  vector  tangent  to  the  yield  circle  at  the  current  force 
point. 


dTn  =  dt  •  n 

(6) 

dTt  =  dt  •  t 

(7) 

dt  =  dTn  n  +  dTt  t  =  dtn  +  dtt 

(8) 

4Ga 

H°  =  ■  is  the  elastic  modulus 

(9) 

K  / 

H  =  H0  (1 - j-jj—  )1/3  is  the  elastoplastic  modulus  corresponding  to  the 

yield  surface  of  radius  K. 

Hp  =  [  1-  (1 - nr-)^3]  *s  tanSent^  elastoplastic 

modulus. 

f  =  coefficient  of  friction  of  the  material  of  the  spheres. 


a  =  (BNR)1/ 3,  with  B  = 
R  =  radius  of  the  spheres. 
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Hardening  Rule: 

In  3— D  force  space,  the  axis  of  any  yield  surface  translates  without  rotation  in 
such  a  way  that  it  remains  always  parallel  to  the  N  axis: 

d t  =  d?  -  fdNn  +  (K  +  fdN)dn  (10) 

where  d£  is  the  translation  of  the  center  of  the  "yield"  surface  on  the  t— plane,  K  is 
the  radius  of  the  yield  surface  on  the  ir— plane  and  dn  is  the  change  of  the  direction 
in  the  unit  vector  n  normal  to  the  surface  on  the  ir~ plane  of  the  current  force  point. 


2.1.1  Verification  of  the  Model 

The  constitutive  relation  briefly  described  above,  implemented  in  computer 
program  CONTACT  (Seridi  and  Dobry  1984,  Dobry  et  al.  1989),  was  subsequently 
tested  to  verify  that  it  reproduces  accurately  the  analytical  solutions  obtained  by 
Mindlin  and  Deresiewicz  (1953).  In  the  numerical  solution,  the  elastic  properties  of 
quartz  were  used  as  input  to  the  program  (Gs  =  32.9  GPa,  i/s  =  0.15  and  f  =  0.5, 
White  1965),  where  Gs  is  the  shear  modulus  and  the  Poisson’s  ratio  of  the 
material  of  the  spheres.  Because  of  the  complexity  of  some  of  the  loading  cases, 
only  six  of  the  most  easily  visualized  cases  will  be  presented  herein.  The 
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force— deformation  curves  from  the  analytical  results  of  Mindlin  and  Deresiewicz 
(1953)  are  reproduced  separately  in  Fig.  10,  while  the  corresponding  numerical 
results  obtained  with  program  CONTACT  are  shown  in  Figs.  11  to  16. 


2.2  Differential  Stress— Strain  Relationships  for  Regular  Arrays  of  Spheres 

Once  the  solution  of  the  contact  problem  is  known,  it  can  be  used  to  develop 
incremental  stress— strain  constitutive  laws  for  regular  arrays  of  spheres.  These 
incremental  stress— strain  relations  can  be  integrated  along  specific  stress  paths  in 
order  to  determine  the  large  strain  response  of  a  regular  packing  of  spheres  under  a 
number  of  loading  conditions. 

Based  on  the  contact  solutions  outlined  in  Section  2.1,  it  is  possible  to  develop 
these  incremental  stress-strain  laws  with  the  help  of  simple  geometric 
considerations  if  the  array  is  statically  determinate,  or  by  using  geometric  and 
compatibility  relations  if  the  array  is  statically  indeterminate.  Such  differential 
stress— strain  relationships  have  been  developed,  among  others,  for  the  simple  cubic 
array  (sc)  shown  in  Fig.  17a  (Deresiewicz  1958),  for  the  body  centered  cubic  array 
(bcc)  of  Fig.  17b  (Petrakis  and  Dobry  1986)  and  for  the  face  center  cubic  array  (fee) 
of  Fig.  17c  (Duffy  and  Mindlin  1957).  These  three  arrays  have  a  number  of  contacts 
per  particle,  or  coordination  number,  CN  of  6,  8,  and  12  respectively.  Other  regular 
arrays  for  which  differential  stress— strain  laws  are  available  include  the  hexagonal 
close-packed  (Duffy  1959)  and  the  cubical-tetrahedral  and  tetragonal  sphenoidal 
arrays  (Makhlouf  and  Stewart  1967). 
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The  general  differential  stress-strain  relationship  for  regular  arrays  of  spheres 
under  isotropic  pressure  is  of  the  form: 


dtTij  —  Sjjitidfki)  or  dejj  —  Cijkido\.i 


(11) 


where  Sijki  is  the  stiffness  and  Cijki  the  compliance  matrix. 

In  what  follows,  the  nonzero  coefficients  of  Sijki  or  Cijki  are  provided  from  the 
studies  listed  above,  for  the  simple  cubic,  body  centered  and  face  centered  cubic 
arrays. 

1.  Simple  Cubic  Array,  (sc,  CN  =  6): 

Snn  =  S2222  =  S3333  =  (|)  /  (1  —  ^s)'2/3  (0oGs)  (12) 

S 1212  =  S13ts  =  S 2323  =  (|)1/3  2[l  Z  -  -  (^oGs2)1/3  (13) 


2.  Body  Centered  Cubic  Array  (bcc,  CN  =  8): 


ani 


=  C2222  =  C3333  = 


(4v^Gs£t0)  ^ 


[1  -  ,s)2/3  + 


2-^s  1 

IZW*1 


(1-^s) 


(14) 
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'1122 


=  C 1 133  —  C2233  — 


75 


1 


(4v/3G  s<70) 


jp  [(2(1  —  ^s)  ^ 


2  -  v  s 
(1  -  vs) 


(15) 


'1212 


=  C 13 13  —  C2323  —  — 


v/3  (4i/3Gsa0) 
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2  —  Vs 
(1  -  "s) 


(16) 


3.  Face  Centered  Cubic  Array  (fee,  CN=12): 


Sim  =  S2222  =  S3333  = 


r  3Gsa0  -I  *■/  ®  4-3 
2(  l-^s)2 


(17) 


Si:22  =  S  H33  —  S2233  — 


r  3Gscr0  i1/3 


(18) 


S 1212  =  s  13 13  =  S2323  = 


r  3Gsi7o  i1/3  4—3 

2(  1  -Vs)2  {2-Vs) 


(19) 


In  Eqs.  12-19,  Gs,  vs  are  the  shear  modulus  and  Poisson’s  ratio  of  the  material  of 
the  spheres,  respectively,  and  cr0  is  the  (macroscopic)  isotropic  pressure  applied  to 
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the  regular  array.  It  should  be  noted  that  for  a  given  a0  Eqs.  11—19  describe 
linearly  elastic,  anisotropic  continua.  This  is  a  result  of  the  crystal  anisotropy 
inherent  to  any  regular  array  structure.  Uniform  sand,  however,  is  more  or  less 
isotropic  under  isotropic  pressure  (Knox  et  al.  1982,  Kopperman  et  al.  1982,  Lade 
and  Nelson  1987),  and  since  the  objective  herein  is  to  model  the  behavior  of  granular 
soil,  it  is  appropriate  to  define  the  conditions  for  which  the  stiffness  matrices  of  the 
arrays  become  isotropic.  From  inspection  of  Eqs.  12—19,  it  can  be  seen  that  the 
necessary  and  sufficient  condition  for  these  three  media  to  be  isotropic  under 
isotropic  pressure  is  that  us  =  0  for  the  material  of  the  spheres.  (Duffy  1959, 
Petrakis  and  Dobry  1986).  The  elastic  constants  of  quartz  are  Gs  2  32.9  GPa  (4780 
ksi)  and  vs  =  0.15  (White  1965).  Moreover,  it  has  been  shown  (Duffy  1959)  that 
even  when  the  degree  of  anisotropy  is  largest  ( t/a  =  0.5)  the  difference  between  the 
moduli  of  the  array  in  different  directions  is  not  more  than  3.6%  to  9%,  depending 
on  the  array  considered. 

Consequently,  us  =  0  is  a  reasonable  first  approximation  for  quartz  spheres 
forming  regular  cubic  arrays.  If  vs  —  0  is  assumed,  Gmax  =  ^  S 1212  =  ^  S1313  =  ^ 
S2323  =  j  Sun,  independent  of  direction  in  the  three  arrays.  For  given  values  of  Gs 
and  (Tq,  Gmax  is  essentially  proportional  to  the  coordination  number,  CN,  as 
illustrated  by  Fig.  18;  the  same  trend  has  also  been  reported  by  Yanagisawa  (1983). 
This  is  very  interesting,  as  it  shows  that  the  stiffness  of  a  regular  array  is  directly 
controlled  by  the  number  of  load-transmitting  interparticie  contacts,  and  it 
suggests  that  an  increased  number  of  contacts  works  in  a  way  similar  to  the  increase 
in  stiffness  of  a  system  of  springs  in  parallel  when  new  springs  are  added.  Of  course, 
as  CN  increases,  the  regular  array  also  gets  denser  and  its  void  ratio,  e,  decreases, 
and  therefore  Fig.  18  can  also  be  interpreted  as  showing  that  Gmax  increases  as  e 
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decreases.  This  inverse  relation  between  Gmax  and  e  is  especially  useful  in  actual 
sands  and  random  arrays,  where  e  can  be  easily  obtained  while  CN  may  be 
impossible  to  determine. 

Based  on  the  expression  for  Graax  with  i/s=  0,  it  is  possible  to  compute 
velocities  of  shear  and  compressional  waves  propagating  through  the  above  three 
regular  packings,  for  wave  lengths  significantly  larger  than  the  radius  of  the  spheres. 
The  shear  wave  velocity  through  an  isotropically  loaded  regular  array  velocity  is 
Vs  =  \JGmAxj p,  where  p  is  the  mass  density  of  the  medium.  In  the  same  way  as 
Gmax,  Vs  becomes  larger  as  CN  increases  and  e  decreases.  A  plot  of  Vs  in  regular 
cubic  arrays  of  quartz  spheres  versus  void  ratio,  e,  for  a  given  isotropic  pressure,  cx0, 
is  presented  in  Fig.  19a,  while  Fig.  19b  shows  the  experimental  curves  for  sands 
reported  by  Hardin  an  Richart  (1963).  While  the  plots  are  qualitatively  similar,  the 
analytically  computed  wave  velocities  for  the  regular  arrays  are  two  to  three  times 
higher  than  those  in  the  actual  soils. 

It  is  hypothesized  herein  that  the  difference  between  Figs.  19a  and  19b  is 
related  to  the  dependence  of  shear  stiffness  on  coordination  number  shown  in  Fig. 
18.  That  is,  the  dependence  of  Gmax  and  Vs  on  the  void  ratio  in  real  soils  is 
explained  by  the  increase  in  the  number  of  intergranular  contacts  as  void  ratio 
decreases.  This  hypothesis  has  also  been  made  by  other  authors  (Deresiewicz  1958, 
Ko  and  Scott  1967).  Furthermore,  while  the  empirical  relations  in  sand  give  Gmax  = 
Af(e)<ro^2  (or  V8  <x  where  f(e)  is  a  function  of  the  void  ratio  e  and  A  is  a 

1/3 

const&nt,  the  3n3lytiC3l  expressions  plotted  in  Fi^.  193  ^ive  Bf(e)  a0  '  (or 

V8  *  aQ  '  ).  These  higher  values  of  moduli  (and  velocities)  obtained  analytically  as 
compared  to  measurements,  and  the  difference  between  the  1/4  and  1/6  power  in 
the  exponent  of  the  confining  pressure  for  Vs,  are  well  known  and  have  been 
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observed  in  the  past  (Duffy  and  Mindlin  1957,  Deresiewicz  1958).  Figure  20 
presents  a  comparison  between  analytically  obtained  and  experimental  values  of  rod 
wave  velocity,  Vr  =  v/Emax/p,  where  Emax  is  the  small  strain  Young’s  modulus, 
versus  confining  pressure,  a0,  in  rods  constructed  of  equal  steel  spheres  assembled  in 
the  fee  packing.  It  should  be  noted  that  as  the  pressure  increases,  not  only  the 
dependence  of  Vr  on  the  pressure  changes  from  a0  1  to  a0  1  ,  but  also  the 
experimental  results  approach  the  theory. 

Duffy  and  Mindlin  and  Deresiewicz  have  provided  an  explanation  for  the 
results  of  Fig.  20,  which  is  consistent  with  the  argument  advanced  herein  to 
interpret  Figs.  19a  and  19b.  As  the  particles  used  for  the  experiments  of  Fig.  20 
were  unequal  in  size  and  not  perfectly  spherical,  they  had  a  number  of  contacts 
which  transmitted  little  load  (spheres  barely  touching  or  touching  on  asperities)  or 
no  load  (spheres  not  touching).  Consequently,  the  measured  stiffness  of  the  actual 
medium  was  less  than  predicted  and  the  dependence  of  the  velocity  on  the  pressure, 
(T0,  was  not  1/6.  When  cr0  increased,  the  barely  touching  contacts  started 
transmitting  their  full  load  and  new  contacts  appeared.  As  a  result,  the 
experimental  values  of  Vr  approached  the  theoretically  predicted  ones.  This 
increase  in  stiffness  and  change  of  the  dependence  of  the  velocity  on  the  exponent  of 
the  pressure  from  1/4  to  1/6  can  be  explained  by  this  initial  increase  and  subsequent 
stabilization  of  the  number  of  effective  contacts  with  pressure,  including  the 
influence  of  the  roughness  of  the  actual  sphere  surface.  The  above  phenomena  are 
illustrated  by  the  two  sets  of  data  in  Fig.  20,  corresponding  to  high  and  low 
tolerance  steel  spheres.  In  the  low  tolerance  spheres,  having  a  greater  variation  of 
spheres’  sizes,  the  velocity  is  lower  but  it  increases  faster  with  <j0  at  low  pressures 
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than  for  the  high  tolerance  spheres.  Both  low  and  high  tolerance  spheres  approach 
the  theoretical  line  and  the  slope  1/6  at  high  pressures  (Deresiewicz  1958). 

Regular  arrays  spheres  capture  some  aspects  of  the  actual  sand  behavior, 
and  some  analytically  obtained  results  are  in  good  qualitative  agreement  with 
laboratory  measurements  on  granular  soils.  As  shown  previously,  they  can  also  be 
made  to  behave  as  a  linear  isotropic  elastic  solid  once  vs  =  0  is  assumed.  However, 
each  one  of  the  three  regular  arrays  has  a  fixed  value  of  void  ratio,  e  (see  Fig.  18) 
and  thus  cannot  simulate  a  sand  of  arbitrary  density.  In  what  follows,  an  improved 
random  model  of  granular  soil  based  on  combining  these  three  regular  arrays  is 
proposed  and  some  of  the  properties  of  this  regular /random  medium  are  compared 
to  measurements  in  actual  sand. 


2.3  The  Self  Consistent  Method 


One  of  the  most  common  procedures  used  to  describe  the  behavior  of 
multiphase  (composite)  media  is  the  Self  Consistent  Method  (Hill  1965,  Budiansky 
1965). 

This  approach  was  first  applied  by  Hershey  (1954)  and  Kroner  (1958)  to  model 
the  behavior  of  isotropic  and  anisotropic  polycrystalline  media.  Hill  (1965)  and 
Budiansky  (1965)  improved  the  method  and  applied  it  to  the  study  of  multiphase 
media;  they  also  developed  the  basic  procedure  adopted  herein.  The  improved 
Hill-Budiansky  method  yields  an  estimate  of  the  macroscopic  elastic  constants  of  a 
multiphase  medium  which  is  an  aggregate  of  isotropic,  linearly  elastic  materials. 
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The  medium  is  assumed  to  consist  of  continuous,  irregular  regions  containing 
the  constituents,  like  those  of  a  polycrystal.  Furthermore,  the  shape  of  these  regions 
is  assumed  not  to  deviate  much  from  ellipsoidal.  The  spatial  distribution  of  the 
microscopic  components  is  such  that  this  "composite"  medium  can  be  assumed  to  be 
macroscopically  isotropic  and  homogeneous.  The  medium  has  N  constituents  and  a 
total  volume  V.  The  volume  of  the  ith  phase  is  Vi  and  the  volume  concentration  of 
the  phase  is  Ci  =  Vi/V.  It  should  be  pointed  out  that  in  the  limiting  case  of  many 
small  concentrations,  Cj,  C2,  C3,  ...c^  _j,  the  first  N— 1  phases  will  tend  to  appear  as 
individual  inclusions  in  a  matrix  which  in  turn  is  the  Nth  phase. 

In  order  to  obtain  the  macroscopic  elastic  constants  of  the  medium,  K*,  a 
uniform  macroscopically  stress  field  a?j,  is  applied  to  the  medium.  Then,  assuming 
that  the  shape  of  the  "inclusions"  is  elliptical,  and  using  the  solution  for  an  elastic 
isotropic  inclusion  inside  an  elastic,  isotropic  medium  (Eshelby  1957),  the  elastic 
field  is  computed  for  each  of  the  phases.  Finally,  in  order  to  determine  the  effective 
moduli,  G*,  */*,  or  K*,  the  total  strain  energy  calculated  in  terms  of  the  individual 
phase  properties,  Gj,  Kj,  is  equated  to  that  of  the  macroscopic,  yet  unknown, 
medium  with  constants  G*  and  K*. 

For  example,  to  compute  the  macroscopic  shear  modulus,  G*  =  r0/ 7,  where 
r°  =  rijdV  is  the  macroscopically  applied  uniform  shear  stress,  see  Hill  (1965), 

J  V 

and  7  is  the  average  value  of  the  engineering  shear  strain,  7ij,  over  the  total  volume, 
V),  the  strain  energy,  E,  is  evaluated  in  terms  of  the  macroscopic  material 
properties  (Budiansky  1965): 

7«dv  (20) 
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In  terms  of  the  individual  phase  properties,  Eq.  20  can  be  rewritten  as  follows 
(Budiansky  1965): 


E  = 


y-o  t 


a 


N 


dV  + 


ij 


_ll 


N 


)dV 


(21) 


which  in  terms  of  the  N  individual  constituent  properties  yields 


N 


(22) 


and  finally, 


E  = 


N— 1 

+  l 

i  =1  w 


(23) 


Comparing  Eqs.  20  and  23,  we  see  that  the  macroscopic  shear  modulus  is  given  by 


N— 1 

=  l  Ci(1 

i  =1 


(24) 


The  expression  for  the  macroscopic  bulk  modulus,  K*,  can  be  obtained  through  a 
similar  calculation 


N— 1 


1  _  1  , 


I  C<(1 


i  =1 


(25) 
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G*,  K*  are  the  unknown  shear  and  bulk  moduli  of  the  medium,  Gj,  Ki  (i  =  1,  ..n) 
are  the  corresponding  moduli  of  phase  i,  Ci  =  Vi/V  is  the  volume  concentration  of 
phase  i,  T,  7vi  the  shear  and  volumetric  strains  of  this  ith  constituent,  and  r°  and 
a°  are  the  applied  shear  and  hydrostatic  stresses  at  the  boundary  of  the  medium. 

Eshelby  (1957)  found  that,  in  the  case  of  an  ellipsoidal  inhomogeneity,  the 
elastic  field  within  the  phase  is  uniform  and  the  corresponding  uniform  strains  are 
given  by: 


7i  = 


_ _ 

G*+/J*(Gi-G*) 


(26) 


0 

tvi  =  K*+a*{Ki-K*) 


(27) 


After  substituting  Eqs.  26  and  27  into  24  and  25,  and  then  by  "smearing  out",  that 
is  substituting  the  properties  of  the  matrix  which  surrounds  every  "inclusion" 
(constituent)  with  those  of  the  resultant,  yet  unknown,  medium,  Eqs.  24  and  25 
simplify  and  become: 


N 

l  _ 

i+1  1  +  P*  (£tt  —  1) 


Ci 

ST 


=  1 


(28) 


N 

i=i 1  +  a*  (r 


C  i 

K7 


(29) 


where  a*  and  b*  are  constants  depending  on  the  shape  of  the  inclusion  (Eshelby 
1957).  For  spherical  inclusions: 


a 


*  _ 


l+i/* 

3(1”"*)' 


4—5  U*) 


and  the  Poisson’s  Ratio  of  the  medium,  u*,  is 


_  K*  -  2G* 

^  -  6K*~+~  2G* 


(30) 

(31) 


(32) 


Finally,  the  elastic  constants  of  the  macroscopic  medium,  G*,  K*,  are 
obtained  by  sol  ring  Eqs.  28—32  simultaneously.  The  solution  obtained  through  this 
method  always  lies  between  the  Voigt  and  Reuss  bounds  (spatial  averages  of  the 
Stiffnesses  and  Compliancies,  respectively)  and,  under  certain  conditions,  between 
the  more  strict  Hashin  and  Shtrickmann  (1963)  bounds  (Hill  1965,  Budiu.nsky  1965). 


2.4  A  Model  of  a  Random  Array  of  Eoual  Spheres 

Smith  et  al.  (1929)  observed  that  after  shaking  and  tamping  had  been  applied 
to  a  random  arrangement  of  equal  spheres,  the  medium  appeared  to  be  composed  of 
clusters  of  dense  and  loose  regular  arrays.  These  measurements  also  showed  that 
CN  ranged  between  6  and  12  contacts  per  sphere,  which  corresponds  to  the 
theoretical  range  for  regular  arrays. 
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Additional  experimental  work  by  Bernal  and  Mason  (1960),  Bernal  et  al. 
(1964),  Scott  (1960),  Davis  and  Deresiewicz  (1974),  Shahinpoor  and  Shahrpass 
(1982),  and  Finney  (1985),  has  confirmed  that  both  2— D  and  3— D  random 
assemblages  of  equal  spheres  tend  to  crystallize.  Consequently,  it  is  generally 
accepted  that  an  assemblage  of  identical  spheres  can  be  modelled  by  a  combination 
of  regular  arrays,  as  idealized  in  Fig.  21a,  where  each  polyhedron  contains  a  regular 
packing.  The  macroscopically  observed  porosity  would  then  be  obtained  by  the 
appropriate  combination  of  the  porosities  of  the  packings  weighed  by  their  volume 
fractions,  Ci  (Finney  1983,  Backman  et  al.  1983). 

Therefore,  the  authors  developed  a  model  of  a  random  array  of  equal  spheres 
consisting  of  randomly  distributed  clusters  of  N  regular  packings  such  as  the  three 
cubic  arrays  discussed  previously.  Furthermore,  every  regular  cluster  is  assumed  to 
be  isotropic  (that  is,  for  cubic  arrays,  the  Poisson’s  ratio  of  the  spheres,  ua  =  0). 
Each  of  these  clusters  is  considered  to  be  an  inclusion  in  the  macroscopic  medium, 
such  as  shown  by  the  observation  of  an  actual  random  2— dimensional  array  of 
spheres  in  Fig.  21b.  where  regions  of  loose  and  dense  packings  can  be  clearly  seen. 
One  could  imagine  the  2D  medium  of  Fig.  21b  to  be  a  cross  section  of  the  3D 
medium  of  Fig.  21a,  which  would  be  composed  of  an  appropriate  combination  of  the 
three  regular  arrays:  fee  (dense  clusters)  bcc  (less  dense  clusters)  and  sc  (loose 
clusters).  In  this  case,  N  =  3. 

As  a  second  step,  a  suitable  combination  of  volume  concentrations  Ci  (i  =  1, 
...N),  which  yields  the  desired  porosity  of  the  random  array,  needs  to  be  defined  for 
N  regular  arrays.  For  this,  the  model  proposed  by  Shahinpoor  (1981)  is  used.  In  his 
work,  Shahinpoor  derived  an  analytical  expression  for  the  spatial  probability  density 
function,  pdf,  of  the  void  ratio,  p(e),  in  a  three-dimensional  random  packing  of 
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equal  steel  spheres,  which  he  later  confirmed  experimentally  in  two  dimensions  by 
means  of  an  optical  scanning  technique  (Shahinpoor  and  Shahrpass  1982). 

Finally,  once  the  volume  fractions,  Ci,  have  been  determined  for  the  regular 
arrays,  the  Self  Consistent  method  is  used  to  average  the  properties  of  the  N  regular 
arrays  within  the  (unknown)  macroscopic  medium,  which  is  compressed  at  its 
boundary  by  a  hydrostatic  pressure,  a°0,  and  subsequently  loaded  by  a  small 
increment  of  stress,  &o\y  Of  course,  for  these  granular  arrays  the  Self-Consistent 
method  must  be  modified  to  account  for  the  effect  of  a0  on  the  elastic  properties  of 
both  the  medium  and  regular  array  clusters. 

The  first  step  in  the  application  of  the  proposed  method  to  random  packings  of 
spheres  is  to  define  the  phases  or  constituents.  Here,  the  sand  will  be  assumed  to 
consist  of  the  three  regular  arrays  of  spheres:  sc,  bcc,  and  fee.  The  expression  of  the 
pdf,  p(e),  is  (Shahinpoor  1981): 


where  A  is  a  constant  defined  by  the  following  equation 


and  enim  emax  are  the  minimum  and  maximum  values  of  the  void  ratios  of  the 
constituents  (efflin  =  0.32,  eraax  =  0.91  for  identical  spheres)  and  e  is  the  mean  value 
of  the  distribution. 
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Since  the  mean  value  of  the  pdf  of  the  porosity,  p(n),  n,  coincides  with  the 
macroscopic  average  (measured  value)  of  the  porosity  of  the  granular  medium, 
whereas  the  mean  void  ratio  does  not,  the  pdf  of  the  void  ratio,  p(e),  is  transformed 
into  p(n)  by  using  the  transformation  p(n)  =  lde/dnjp(e)  (Benjamin  and  Cornell 
1970).  Then  the  appropriate  p(n)  is  chosen  so  that  n  coincides  with  the 
macroscopically  measured  value.  The  continuous  function  p(n)  is  subsequently 
lumped  at  the  three  porosities  of  the  regular  arrays  and  thus  the  volume 
concentrations  are  determined.  If  now  it  is  assumed  that  the  above  three  regular 
packings  form  clusters  ("inclusions")  within  the  macroscopic  medium,  which  are 
approximately  ellipsoidal  in  shape  (spherical,  in  this  instance  for  simplicity)  the  Self 
Consistent  Method  described  by  Eqs.  24—32  can  be  applied.  This  case  of  a  granular 
medium  is  rather  complex  since  the  properties  of  the  constituents,  Gi,  Ki,  as  well  as 
of  the  whole  medium,  G*,  K*,  are  pressure  dependent.  Furthermore,  the  pressure 
taken  by  each  inclusion,  al,  which  influences  directly  its  elastic  constants  (see  Eqs. 
11—18),  depends  on  the  moduli  of  both  the  inclusion  and  the  macroscopic,  yet 
unknown,  medium.  Once  again,  the  pressure  aQ  inside  inclusion  i  can  be  computed 
from  the  Eshelby  solution  as  follows: 


i 

(Jo  = 


0  Ki  r 
°0  [■ 


1  +  a* 
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(35) 


Finally,  the  problem  reduces  to  the  simultaneous  solution  of  Eqs.  28—32  and  35  in 

which  the  values  of  Gi  =  Gi(ao)  and  Ki  =  Ki(cro)  are  given  by  Eqs.  11-18  and  the 

« 

volume  concentrations,  Ci,  from  the  lumped  p(n)  for  i  =  1,  2  and  3  corresponding  to 
the  sc,  bee  and  fee  arrays  respectively. 
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2.5  Application  of  the  Self— Consistent  Method  to  Uniform  Rounded  Sand 

The  proposed  model  was  evaluated  by  predicting  the  small  strain  shear 
moduli,  Gmax,  of  a  rounded,  uniform,  quartz  sand  which  had  been  subjected  to 
isotropic  pressure  in  the  resonant  column  device  and  was  then  heavily  precycled 
(Fig.  2).  This  was  achieved  by  using  as  input  the  elastic  properties  of  quartz  for  the 
individual  spheres  (grains):  Es  =  75.8  GPa  (from  Gs  =  32.9  GPa  and  us  —  0.15, 
White  1965).  However,  in  order  to  ensure  that  the  individual  arrays  as  well  as  the 
macroscopic  medium  are  isotropic  under  isotropic  pressure,  a  value  of  us  =  0  was 
assumed  here  for  the  Poisson’s  ratio  of  the  spheres  and  used  with  Es  =  75.8  GPa. 
The  actual  value  of  vs  is  very  small  already,  and  as  mentioned  previously,  this 
substitution  does  not  affect  the  values  of  the  moduli  of  the  arrays  in  different 
directions  by  more  than  about  2-3%. 

The  computed  values  of  Gmax  =  G*  appear  in  Fig.  22  (continuous  lines) 
plotted  versus  the  isotropic  confining  pressure  acting  on  the  medium,  <Xq,  for  e  = 
n/(l-n),  equal  to  0.46,  0.54  and  0.58.  The  corresponding  values  of  bulk  modulus, 
K+  =  tfo/fy,  were  ^s0  computed,  and  Fig.  23  contains  a  plot  of  the  confining 
pressure,  a°0)  versus  the  resulting  average  volumetric  strain,  ev,  predicted  by  the 
model  for  e  =  0.54.  The  volumetric  strain  was  derived  from  this  bulk  modulus,  K*, 
for  different  0%:  lv  =  al/K*. 

Figure  22  also  contains  experimental  data  from  precycled  Ottawa  C— 109  sand 
where  millions  of  cycles  were  applied  to  isotropically  consolidated  hollow  cylindrical 
specimens.  This  is  the  same  research  previously  discussed  in  connection  with  Fig.  2 
(Drnevich  1967).  The  data  in  Fig.  22  are  for  the  same  void  ratios  used  in  the 
analytical  simulations,  with  a  number  of  cycles,  N,  between  1  X  108  and  22  X  10®, 
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and  an  applied  cyclic  torsional  shear  strain,  7^,  between  3  X  10'4  and  6  X  10*4. 
This  precycling  increased  the  sand  stiffness  more  than  three  times,  while  at  the 
same  time  the  void  ratio  stayed  about  constant  and  equal  to  the  initial,  virgin 
specimen  (dotted  line  in  the  figure).  Furthermore,  the  precycling  modified  the 
dependence  of  Gmax  on  pressure,  a0^  \  from  about  to  (Fig.  22c),  and 
made  the  experimental  results  approach  closely  the  theoretical  values  after  a  very 
large  number  of  cycles.  Those  trends  provide  strong  confirmation  to  the  hypothesis 
that  sand  behaves  as  predicted  by  the  theory,  once  the  theoretical  number  of  load 
transmitting  contacts  is  reached,  in  this  case  as  a  result  of  many  cycles  of  shear 
precycling. 

In  Fig.  23,  an  additional  hydrostatic  compression  was  applied  to  a  specimen 
after  precycling  it  1  X  108  times  at  a  cyclic  shear  strain  of  6  X  10  *4  (square  data 
points).  Here,  the  stiffness  of  the  medium  approached  the  theoretical  values  faster 
than  did  the  shear  stiffness  in  Fig.  22,  for  the  same  number  of  cycles.  This 
happened  probably  because,  in  addition  to  the  precycling  there  was  a  subsequent 
increase  in  the  confining  pressure,  which  completed  the  formation  of  contacts 
started  by  the  repeated  shearing. 

Additional  experimental  data  from  resonant  column  tests  on  30-40  Ottawa 
C— 109  sand,  this  time  on  solid  specimens,  covering  a  wide  range  of  densities  from 
dense  to  medium  loose,  ar*  presented  in  Fig.  24  (Song  and  Stokoe,  1987).  The 
strain  amplitude  in  these  experiments  was  7^=  6  X  10'4.  It  cam  be  clearly  seen 


(*)  <70  and  <7o  are  both  used  interchangeably  herein  in  the  text  and  figures  to 

denote  the  macroscopic,  isotropic  pressure  applied  to  the  medium. 
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that  while  again  the  value  of  the  shear  modulus,  Gmax,  increases  dramatically  with 
number  of  cycles  (Fig.  24a),  the  void  ratio  remains  essentially  constant  (Fig.  24b). 
It  becomes  almost  certain,  then,  that  this  increase  in  the  stiffness  is  not  caused  by 
any  great  change  in  the  particulate  structure  of  the  medium,  which  would  have 
manifested  itself  through  a  change  in  the  macroscopic  void  ratio,  but  by  the  much 
subtler  change  associated  with  an  increase  of  the  number  of  contacts  at  constant 
void  ratio,  as  one  could  deduce  from  Figs.  20  and  24.  With  the  continuous  cycling 
at  small  strains,  probably  not  only  new  contacts  were  created  between  the  sand 
grains,  but  old  "dead"  ones  became  load-transmitting  ("live"),  resulting  in  the 
excellent  agreement  between  theory  and  experiment  in  Fig.  22  after  many  cycles. 

Finally,  Fig.  25  presents  the  variation  of  shear  modulus  at  small  strains,  Gmax, 

with  number  of  contacts  per  sphere  CN  and  with  cr0  calculated  from:  i)  the  three 

regular  cubic  arrays,  that  is  from  Fig.  18;  ii)  the  Self  Consistent  Method;  and  iii) 

the  analytical  solution  of  Walton  (1987).  Walton,  by  considering  the  pressure 

dependent  normal  and  tangential  compliances  at  the  interparticle  contacts  derived 

*  * 

the  following  expressions  for  the  two  elastic  moduli  (Lame  constants),  A  and  n  = 
Gmax i  of  a  random  packing  of  equal  elastic  spheres  with  an  infinite  coefficient  of 
friction  under  an  isotropic  pressure  a0: 
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where  <^?=(  1 — n),  and  n  is  the  average  porosity  of  the  medium,  B  and  C  are  constants 
which  depend  on  the  material  properties  of  the  spheres,  Gs  and  va,  and  CN  is  the 
average  coordination  number. 

The  final  result  is  the  very  consistent  plot  of  Fig.  25,  with  the  Self  Consistent 
Method  providing  a  prediction  in  excellent  agreement  with  the  analytical  results  by 
Walton  (in  Eqs.  35  and  36,  us  was  set  to  zero  and  the  coordination  number  was 
given  by  this  Self  Consistent  scheme).  Thus,  the  earlier  hypothesis,  extrapolated 
from  the  regular  arrays,  that  the  small  strain  shear  modulus,  Gmax,  of  a  random 
array  of  equal  spheres  is  essentially  a  linear  function  of  the  average  number  of 
contacts  per  particle  is  confirmed  and  justified  by  two  different  analytical 
approaches. 


2.6  A  Two-Dimensional  Numerical  Model  of  Granular  Soil  at  Small  Strains 

In  the  previous  section,  a  closed  form  solution  was  obtained  for  the  elastic 
constants  of  a  random  aggregate  of  equal,  rough,  elastic  spheres  having  an  arbitrary 
macroscopic  void  ratio  and  subjected  to  isotropic  loading.  This  was  done  through 
the  use  of  the  Self-Consistent  Model  and  on  the  assumption  that  the  array  is 
composed  of  different  phases,  and  the  results  contributed  to  an  improved 
understanding  of  the  small  strain  behavior  of  sands  under  isotropic  conditions.  The 
fact  that  this  analytical  solution  was  obtained  with  relatively  small  effort  should  be 
attributed  to  the  high  level  of  symmetry  of  such  a  system  under  isotropic  pressure. 

Since  the  Self  Consistent  Method  was  developed  to  determine  the  elastic 
constants  of  a  multiphase  medium,  it  is  not  suitable  for  the  determination  of 
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elastoplastic  properties.  Although  there  have  been  modifications  proposed  to  adapt 
the  Self  Consistent  Method  to  nonlinear  behavior,  these  adaptations  would  not  work 
in  the  case  of  a  granular  medium,  since  particle  sliding  changes  the  material 
configuration  and  thus  it  also  changes  the  distribution  of  the  various  volume 
fractions.  Moreover,  the  symmetry  of  the  constituents  would  change  dramatically 
with  loading  and  the  problem  would  become  intractable.  To  overcome  this,  it  was 
decided  to  attempt  a  numerical  model,  which  although  containing  assumptions 
similar  to  those  of  the  Self  Consistent  Method,  could  perform  numerically  the 
integration  of  the  equations  to  larger  strain  levels.  This  model  is  similar  to  those 
proposed  by  a  number  of  researchers  (Budiansky  and  Wu  1962,  Lin  and  Ito  1965, 
1966)  as  discussed  in  Section  1.1,  in  that  it  is  composed  of  regular  packings  of 
spheres  randomly  oriented  so  as  to  ensure  a  statistically  isotropic  medium.  The 
properties  of  the  packing  to  be  used  as  basic  "element"  of  this  medium  is  defined  by 
the  contact  law  described  in  Section  2.1.  In  other  words,  the  constitutive  law  of  this 
regular  packing  is  similar  to  those  described  by  Eqs.  12—18,  but  now  the 
compliances  are  given  by  program  CONTACT  described  in  Section  2.1. 


2.6.1  Nonlinear  Finite  Element  Simulations 

In  order  to  model  the  behavior  of  granular  soil  at  small  strains,  a  two- 
dimensional  numerical  model  was  developed  which  calculates  the  response  of  a 
random  aggregate  of  equal,  elastic,  rough  spheres  under  an  arbitrary  boundary  stress 
state,  o jj .  For  this,  a  finite  element  analysis  was  performed  in  which  the  element 
corresponded  to  a  simple  cubic  array  of  equal,  rough  elastic  quartz  spheres.  Each 
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2D  element  contains  an  undetermined  number  of  spheres,  and  is  assumed  to  be 
subjected  to  a  uniform  stress  field  so  that  four  spheres  can  be  used  to  represent  the 
above  element  (Duffy  and  Mindlin  1957).  Figure  26  sketches  this  individual 
element.  The  stress— strain  behavior  of  this  array  under  a  biaxial  state  of  stress,  au 
and  cr22,  followed  by  pure  shear,  appears  in  Fig.  27c.  Figures  27a  and  27b  portray 
the  force— deformation  behavior  of  the  '’weak"  and  "strong"  contacts  for  the  same 
case,  respectively.  Therefore,  the  behavior  of  this  packing  under  biaxial  loading  is 
the  result  of  the  interaction  of  the  two  sets  of  contacts.  In  this  element,  the 
force— deformation  relation  at  each  of  the  contacts  between  spheres  is  given  by  the 
numerical  solution  discussed  in  Section  2.1,  through  program  CONTACT,  and  the 
constitutive  law  is  obtained  by  combining  the  numerical  solutions  through 
geometric  considerations.  The  element  cannot  sustain  a  tensile  stress  normal  to  the 
slip  planes,  since  this  would  imply  that  the  corresponding  contacts  would  cease  to 
exist,  and  as  a  result  the  particles  would  attempt  to  rearrange  and  form  a  new 
packing.  Although  this  rearrangement  obviously  happens  in  actual  sand  aggregates, 
its  simulation  is  too  complex  and  computationally  demanding  to  implement  in  the 
present  finite  element  code,  and  the  decision  was  made  to  make  the  contact  normal 
forces  positive-definite,  that  is  they  can  be  either  positive  or  zero,  but  never 
negative. 

A  subroutine  implementing  the  behavior  of  this  element  was  coded  into  the 
nonlinear  finite  element  program  ABAQUS  (1982)  as  a  user  defined  material 
subroutine  (UMAT).  Finally,  an  incrementally  linear  analysis  was  done  using  eight 
noded  elements  with  reduced  integration. 

Several  media  were  used  in  these  simulations  consisting  of  a  number  of  simple 
cubic  arrays  in  two  dimensions  (a  monolayer  of  equal  spheres  assembled  in  simple 
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cubic  patterns),  oriented  in  a  such  a  way  as  to  resemble  a  statistically  isotropic 
random  aggregate.  As  discussed  above,  each  of  the  elements  is  a  simple  cubic  array, 
randomly  oriented,  and  the  number  of  elements  used  in  each  of  these  media  varies 
from  16  (4  x  4)  to  64  (8  x  8)  as  illustrated  in  Fig.  28.  Therefore,  the  medium  has 
the  same  properties  of  the  simple  cubic  array  in  terms  of  void  ratio  and  coordination 
number.  Although  other  packings  could  in  principle  have  been  used,  it  would  have 
been  more  difficult  (Petrakis  and  Dobry  1987).  It  was  felt  that  the  inability  of 
varying  the  void  ratio  or  the  coordination  number  (number  of  contacts  per  particle) 
did  not  unduly  constrain  the  planned  FE  simulations,  which  were  expected  to 
provide  useful  insight  into  the  behavior  of  a  sand  at  small  strains. 


2.6.2  Monotonic  Finite  Element  Loading  Simulations 

As  a  first  step,  the  isotropy  of  the  media  used  was  verified  by  loading  them 

isotropically  up  to  o°  =  100  KPa,  followed  by  pure  shear  which  was  applied 

o 

incrementally.  This  was  accomplished  by  imposing  on  each  medium  a 
predetermined  direction  of  the  major  principal  stress.  The  values  used  for  the  angle, 
a,  between  the  major  principal  stress  and  the  vertical  direction  of  the  medium  (Fig. 
29a)  were  0°  (compression  in  the  vertical  direction);  22.5°;  45°  (pure  shear  in  the 
vertical  and  horizontal  planes);  67.5°;  and  90°  (extension  in  the  vertical  direction). 
The  results  of  these  simulations  are  shown  in  Fig.  30  as  plots  of  the  applied  deviator 
stress,  a  —  cr°,  versus  the  resultant  shear  strain,  7  =  e*  — e2-  Figure  30  includes 
results  of  all  four  media  composed  of  16  and  64  elements  of  Fig.  28. 
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It  can  be  seen  in  Fig.  30  that  the  medium  is  indeed  isotropic  under  isotropic 
pressure,  as  expected.  Since  Fig.  30  shows  that  the  difference  between  the 
stress-strain  behavior  of  the  16— element  and  64-element  media  is  not  significant,  it 
was  decided  that  for  subsequent  parametric  studies,  as  well  as  for  monitoring  the 
stress-strain  behavior  of  each  element,  any  of  the  less  costly  16— element  media 
could  be  used  as  representative  of  the  aggregate.  In  another  simulation  using 

Medium  2  with  16  elements,  a  monotonically  increasing  hydrostatic  pressure  was 

o 

applied  up  to  a ^  =  500  KPa,  and  it  was  observed  that  both  the  macroscopic  and  the 
microscopic  (element)  response  exhibits  a  locking  nonlinear  elastic  behavior,  similar 
to  that  observed  in  sand  (Fig.  31,  see  also  Fig.  23). 

Therefore,  the  media  being  used  in  the  numerical  experiments  simulate  well 
important  aspects  of  the  behavior  of  actual,  uniform,  rounded  sand,  by  being 
isotropic  under  isotropic  pressure,  yielding  in  shear  and  locking  under  hydrostatic 
compression.  That  is,  these  "random"  media  have  the  desirable  properties  of  the 
regular  simple  cubic  array  without  its  problematic  aspects. 

In  the  numerical  simulations  of  pure  shear  summarized  in  Fig.  30,  it  was 
observed  that  the  yielding/failure  process  of  the  medium  occurs  in  two  successive 
stages.  In  the  first  stage,  a  growing  number  of  "soft"  elements,  oriented  more  or 
less  parallel  to  the  directions  of  the  applied  shear  stress,  slide  and  this  sliding 
accounts  for  the  increasing  nonlinearity  of  the  curves  in  Fig.  30,  as  the  shear  strain, 
7  =  — 12,  increases  from  0  to  values  around  0.1  X  10*3.  At  these  larger  values  of 

shear  strain,  typically  around  20%  of  the  elements  have  already  slid  (failed).  In  the 
second  stage,  occurring  at  about  7  =  —  e2  =  0.16  x  10'3,  one  or  several  of  the 
"stiff1  elements,  oriented  more  or  less  perpendicular  to  the  direction  of  the  shear 
stress,  and  which  had  not  yet  slid,  tend  to  separate  as  the  normal  force  at  the 
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contact  becomes  zero,  the  corresponding  ratio  shear/normal  force  at  the  contact 
reaches  f,  and  the  element  slides.  This,  of  course,  is  related  to  the  fact  that  the 
normal  contact  force  is  allowed  to  be  zero  but  not  negative.  Once  some  of  the 
"stiff1  elements  fail  due  to  this  tendency  to  separate,  a  growing  number  of  both 
"soft"  and  "stiff"  elements  slide  in  the  next  increment(s)  by  a  combination  of  shear 
stress  increase  and  separation  tendency,  thus  precipitating  the  failure  of  the 
medium,  occurring  at  rf  =  44.8  KPa.  This  value  of  Tf  is  close  to,  but  slightly  less 
than  the  yield  stress  of  the  simple  cubic  array  subjected  to  pure  shear  on  the  sliding 
planes  of  the  array:  r  =  (0.5)(100)  =  50.0  KPa.  This  "failure"  of  the  aggregate, 
defined  here  by  the  sequence  of  phenomena  just  described,  which  at  the  end  result  in 
the  global  stiffness  matrix  of  the  medium  becoming  singular,  is  associated  with  a 
generalized  tendency  of  the  particles  to  slide,  separate  and  rearrange  themselves  into 
more  stable  positions.  This  corresponds  roughly  to  the  changes  in  geometry 
occurring  in  actual  sands  at  the  threshold  strain,  jt  =  0.1  to  0.2  x  10 "3  (Dobry  et  al. 
1982),  as  verified  by  the  fact  that  "failure"  of  the  medium  in  Fig.  6  occurs  at  a  shear 
strain,  7  =  et  — e2  =  0.16  x  10'3. 


2.6.3  Cyclic  Finite  Element  Loading  Simulations 


The  (16-element)  Medium  1  was  also  subjected  to  an  isotropic  stress  of 

a  =300  KPa  followed  by  a  complete  isotropic  cycle  with  an  amplitude  200  KPa, 
0 

that  is  a  was  first  increased  to  500  KPa  and  then  decreased  to  100  KPa,  and  then 

0 

back  to  500  KPa.  The  stress-strain  behavior  calculated  for  the  medium  is 
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pressented  in  Fig.  32,  and  it  shows  a  nonlinear  Mastic  behavior  similar  to  that  of 
actual  sands  (Ko  and  Scott  1967). 

o 

The  same  16-element  medium,  consolidated  isotropically  to  <r  =  100  KPa, 

o 

was  also  cycled  under  pure  shear  conditions,  with  an  amplitude  of  shear  stress,  rc  = 
20,  30,  35,  40  and  43  KPa.  The  hysteresis  loops  for  rc  =  40  and  rc  =  43  KPa 
appear  in  Figs.  33a  and  33b,  respectively. 


2.6.4.  Computation  of  Dynamic  Properties  and  Wave  Propagation  Velocities. 

The  secant  shear  moduli,  G,  obtained  from  the  pure  shear  simulations  at 
different  angles  a  for  the  constant  mean  stress  case  (Fig.  30),  were  normalized  with 
respect  to  the  shear  moduli  at  very  small  strains,  G,^,  obtained  in  the  same 
simulations,  and  the  corresponding  values  of  G/Gm&x  versus  shear  strain  are  plotted 
in  Fig.  34a,  where  they  are  compared  with  the  bounds  proposed  by  Seed  and  Idriss 
(1970)  from  actual  tests  on  sands.  The  corresponding  damping  ratio,  f,  obtained 
from  the  loops  under  cyclic  pure  shear  (Fig.  33)  is  plotted  against  cyclic  shear  strain 
in  Fig.  34b.  Figure  34  corresponds  to  simulations  done  at  a  mean  stress  of  100  KPa. 
Figure  34  also  includes  the  corresponding  curves  for  G/ G^,.  and  damping  ratio  for 
the  case  of  pure  shear  along  the  sliding  planes  of  a  simple  cubic  array,  computed  by 
Dobry  et  al.  (1982).  The  slopes  of  both  the  G/G^  and  damping  ratio  versus  7 
curves  in  Fig.  34  for  the  aggregate  are  flatter  than  that  of  the  simple  cubic  array. 
This  "stiffer"  behavior  of  the  aggregate,  as  compared  to  that  of  a  single  cubic  array, 
should  be  attributed  to  the  interaction  between  the  elements  of  the  medium. 
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Finally,  the  small  strain  constrained  moduli,  D^,  were  computed  in  both 

principal  (vertical  and  horizontal)  directions,  as  follows:  The  medium  was  loaded  to 

o  o 

the  desired  biaxial  stress  ratio,  K  =  a  a  ,  and  then  very  small  stress  increments 

r  2 

o  o 

with  the  appropriate  sign  were  applied,  Act  ,  A  a ^  the  corresponding  differences  in 

strain  were  computed  and,  finally,  the  small  strain  constrained  moduli  of  the 

medium,  Du,  D22  were  calculated  in  both  directions.  The  results  of  this  simulation 

are  shown  in  Figs.  35a  and  35b  as  plots  of  the  normalized  constrained  moduli  of  the 

medium,  D^/D^  and  D^/D^  versus  the  stress  ratio  K  =  a  /a  ,  where  is 
22  22  1  1  11  1  2  ii 

the  constrained  modulus  at  a  given  K,  and  the  corresponding  constrained 

ii 

0 

modulus  under  the  initial  isotropic  stress,  a  .  The  same  figure  includes  also  data 

o 

points  from  a  number  of  measurements  on  sand  in  the  large  cubic  facility  at  the 
University  of  Texas  at  Austin  (Kopperman  et  al.  1982),  which  were  performed 
during  a  test  similar  to  that  defined  for  these  numerical  simulations.  The 
agreement  between  normalized  experimental  results  and  numerical  simulations  in 
Fig.  35a  and  35b  is  excellent.  Consequently,  the  main  conclusion  obtained  from  the 
University  of  Texas  laboratory  results,  that  the  P— wave  velocity  propagating  along 
a  principal  stress  direction  is  only  a  function  of  the  value  of  that  principal  stress  is 
fully  predicted  by  the  numerical  experiments.  Therefore,  as  previously  hypothesized 
by  the  authors,  this  is  an  effect  due  to  the  particulate  nature  of  the  soil,  which  can 
be  explained  and  reproduced  analytically  once  this  particulate  nature  is  taken  into 
account. 
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2.7.  Two-Dimensional  Nonlinear  Distinct  Element  Simulations 


After  the  discussion  on  the  limitations  of  the  nonlinear  finite  element  model 
described  above,  it  was  decided  to  use  a  different  type  of  numerical  simulation, 
capable  of  modelling  granular  soil  subjected  to  large  strain  anisotropic  loading,  in 
which  contact  separation,  sliding  and  rearrangement  of  the  particles  are  allowed. 

There  are  several  approaches  to  model  the  behavior  under  load  of  random 
arrays  of  equal  and  unequal  spheres.  One  of  the  most  widely  used  procedures  is  the 
"Distinct  Element  Method"  (Cundall  and  Strack  1979),  who  have  used  an  explicit 
finite  difference  scheme  to  determine  the  static  response  of  an  array  to  applied 
strains  (program  TRUBAL)  or  to  applied  boundary  displacements  (program 
BALL). 

Program  TRUBAL  (Strack  and  Cundall  1984)  was  modified  at  RPI  by  Ng  and 
Dobry  (1989)  within  another,  NSF  sponsored  project,  by  replacing  the  existing 
arbitrary  linearly  elastic— perfectly  plastic,  non  pressure  dependent,  force 
displacement  relation  at  the  intergranular  contacts,  by  the  more  realistic  and 
rigorous  Hertz— Mindlin  contact  law.  This  was  accomplished  once  more  by 
attaching  program  CONTACT,  described  in  Section  2.1,  as  a  subroutine  in 
TRUBAL  (Table  2).  The  new  two-dimensional  program  was  named  CONBAL-2  ; 
it  is  very  similar  to  three-dimensional  program  TRUBAL,  but  it  has  one  less 
dimension,  due  to  the  increased  computation  and  storage  requirements  imposed  on 
the  program  by  the  nonlinear  subroutine  CONTACT.  Moreover,  program 
CONBAL— 2  has  been  vectorized  to  run  more  efficiently  on  a  supercomputer  with 
vector  facilities.  Another  difference  is  the  absence  of  particle  rotation  in 
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CONBAL— 2,  as  this  causes  as  yet  unsolved  problems  in  the  numerical  simulations. 
This  absence  of  rotation  as  well  as  the  2— D  nature  of  the  simulation  make  the 
random  arrays  of  spheres  somewhat  stiffer  than  actual  3— D  arrays  and  soils;  Ng  and 
Dobry  (1989)  have  accounted  for  this  by  reducing  the  interparticle  angle  of  friction 
from  0.5  to  0.35.  Therefore,  CONBAL— 2  is  very  similar  to  TRUBAL,  except  for 
the  differences  noted  above. 

The  results  of  extensive  parametric  studies  performed  to  check  the  accuracy  of 
this  new  program  appear  in  Ng  (1989).  These  parametric  studies  established  the 
appropriate  range  of  the  necessary  input  variables  (time  step,  strain  rate  and 
damping)  and  confirmed  the  excellent  agreement  between  the  results  of  CONBAL— 2 
and  existing  analytical  results.  At  this  point,  a  copy  of  CONBAL— 2  became 
available  to  the  AFOSR  project  for  research  on  the  stress— strain  behavior  of 
granular  soil. 

As  a  first  step,  it  was  decided  to  start  simulating  phenomena  in  the  small 
strain  range,  and  at  a  later  stage  to  focus  on  the  large  strain,  fully  nonlinear 
inelastic  behavior.  Furthermore,  it  was  also  decided  to  compare  the  results  of 
CONBAL— 2  with  answers  obtained  using  the  analytical  and  numerical  procedures 
presented  in  Section  2,  as  well  as  with  experimental  data  in  the  literature. 

Since  the  interparticle  force— deformation  relationships  developed  by  Cattaneo 
(1938),  Mindlin  (1949)  and  Mindlin  and  Deresiewicz  (1953)  apply  only  to  spheres  of 
the  same  sizes  and  same  properties,  only  random  arrays  of  equal  and  moderately 
unequal  spheres  were  developed  and  studied.  Program  CONTACT,  and  thus 
CONBAL— 2  can,  in  principle  handle  only  random  arrays  of  equal  spheres;  however, 
the  same  program,  with  certain  assumptions,  can  handle  in  first  approximation 
moderately  unequal  spheres  of  the  same  material. 
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2.7.1  Small  Strain  Moduli  of  Isotropically  Compressed  Random  Arrays  of 
Equal  Spheres. 

Two  different  2— D  random  arrays  of  477  identical  particles  having  the 
properties  of  quartz  were  generated  using  CONBAL— 2.  The  first  array  was  very 
loose  with  an  average  coordination  number,  CN,  of  2.1,  while  the  second  was 
medium  dense  with  CN=3.0.  These  arrays  were  then  subjected  to  three  different 
values  of  isotropic  pressure  (<7o=91,  334  and  698  KPa)  without  significant  change  in 
their  average  number  of  contacts. 

The  configuration  of  the  array  with  CN=3.0  subjected  to  cr^  =  =  aQ  =  91 

KPa,  is  shown  in  Fig.  36,  where  the  circles  represent  the  spheres  and  the  rectangles 
the  relative  magnitudes  and  directions  of  the  contact  forces.  There  are  four 
different  rectangle  widths,  with  each  one  of  them  corresponding  to  a  range  of  forces 
between  four  equal  fractions  of  the  maximum  computed  contact  force.  For  example, 
if  the  maximum  contact  force  is  F  KN,  the  narrowest  rectangle  stands  for  the  range 
of  forces  between  0  and  F/4  KN,  the  next  wider  rectangle  for  the  range  of  forces 
between  F/4  and  F/2  KN,  etc.  This  notation  will  be  used  throughout  this  report 
whenever  results  from  the  distinct  element  method  are  presented. 

Figure  37  portrays  some  of  the  micromechanical  statistics  which  are  calculated 
as  a  result  of  the  isotropic  compression  simulation.  Figure  37  includes  the  frequency 
distributions  of  contact  angles,  mobilized  angle  (angle  between  contact  force  and 
contact  normal),  coordination  number  and  contact  force.  As  can  be  seen  in  the 
above  distributions,  the  array  is  not  perfectly  isotropic.  This  is  a  result  of  the  small 
number  of  spheres  used  in  the  periodic  cell,  coupled  with  the  crystalization  which 
occurred  due  to  the  identical  sizes  of  all  477  spheres.  The  combined  effect  of  these 
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two  features  has  increased  the  "characteristic  size"  of  the  smallest  constituent  to 
that  of  a  crystalized  region;  as  a  result  the  array  is  not  statistically  isotropic. 
Another  important  observation  is  the  transmission  of  the  applied  forces  through 
columns  of  particles,  in  which  all  "columns"  experience  the  same  contact  force;  this 
has  caused  the  contacts  to  form  along  preferential  directions  rather  than  uniformly 
in  all  directions. 

These  deviations  from  isotropic  behavior  are  not  that  significant,  however,  and 
the  477— particle  arrays  can  be  considered  to  be  isotropic,  ai  least  in  first 
approximation. 

These  two  arrays  of  477  equal  particles  were  subsequently  used  to  compute  the 
small  strain  shear  modulus  Gmax  under  three  levels  of  isotropic  pressure:  aQ  =  91, 
334  and  698  KPa.  The  shear  modulus  was  computed  by  applying  an  increment  of 
macroscopic  shear  strain,  A7  =  1  X  10'8,  and  computing  Gmax  as  Gmax  =  At  /  A 7. 
The  results  of  this  simulation  appear  in  Fig.  38  (diamonds),  as  a  plot  of  Gmax  versus 
the  Coordination  Number.  In  this  figure  they  are  also  compared  with  the  values 
obtained  from  the  regular  arrays  of  spheres,  the  Self  Consistent  method  and  the 
analytical  expressions  of  Walton  (1987),  already  discussed  in  Section  2.5  and  plotted 
in  Fig.  25.  In  order  to  compare  the  2-D  results  obtained  by  the  Distinct  Element 
Method  with  those  obtained  by  other  3-D  approaches,  it  was  decided  to  adjust  the 
coordination  number  to  account  for  the  third  dimension.  The  only  regular  packing 
of  identical  spheres  where  a  clear  relation  can  be  established  between  the 
coordination  number  and  Gmax  iu  two  and  three  dimensions  is  the  simple  cubic 
array,  which  has  the  same  shear  moduli  in  2  and  3-D.  In  this  regular  array,  CN=4 
for  2— D  and  CN=6  for  3— D.  Therefore  the  values  of  the  coordination  number  in  the 
Distinct  Element  Simulations  of  Fig.  38  were  multiplied  by  6/4=  1.5.  The  final 
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result  is  the  very  consistent  plot  of  Fig.  38  using  points  from  four  different  methods. 
The  plot  essentially  confirms  the  hypothesis  presented  by  Yanagisawa  (1983)  and 
Petrakis  and  Dobry  (1986)  that  the  small  strain  shear  modulus,  Gmax,  of  a  random 
array  of  equal  spheres  is  in  first  approximation  a  linear  function  of  the  average 
number  of  contacts  per  particle. 


2.7.2  Small  Strain  Moduli  of  Anisotropicallv  Compressed  Random  Arrays 

The  same  2— D  array  already  discussed  with  477  equal  elastic,  rough  spheres 
and  coordination  number,  CN=3,  consolidated  under  a0  =  9lKPa  (Fig.  36),  was 
further  loaded  in  program  CONBAL— 2  under  biaxial  compression  to  cm  =  233 
KPa,  while  keeping  constant.  This  was  done  to:  i)  investigate  the  influence  of 
the  magnitude  and  direction  of  the  principal  stress  on  the  constrained  moduli,  D, 
and  also  on  the  velocities  of  P— waves  propagating  through  the  medium,  ii)  interpret 
the  experimental  findings  of  Kopperman  et  al.  (1982),  and  iii)  verify  the  results  of 
the  nonlmear  finite  element  model  presented  in  Section  2.6. 

The  constrained  moduli,  Du  =  p  Vp  were  calculated  in  the  same  way  as  in  the 
nonlinear  finite  element  model  of  Section  2.6,  as  follows:  once  the  desired  stress 
ratio,  K  =  a  Jcj  =  cm/ an,  was  reached  at  specific  points  during  loading,  very  small 
increments,  deu,  de^,  equal  to  1.062  X  10 "6  were  applied  to  the  array  with  the 
appropriate  signs.  The  corresponding  stress  increments,  dan,  da22,  were  computed, 
and  finally,  the  small  strain  constrained  moduli,  =  dau/deii,  were  calculated 
both  in  the  direction  of  the  major  (a  1  =  (722)  and  minor  (<72  =  an)  principal  stresses. 
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The  same  simulation  was  performed  on  a  second  random  array  of  531 
moderately  unequal  quartz  spheres,  having  a  ratio  of  radii,  R1/R2  =  1.5.  In  this 
array  168  spheres  had  a  radius  Rj  and  363  had  a  radius  of  R2.  The  relation  of  sizes 
was  set  to  1.5  so  that  the  Mindlin— Deresiewicz  theory  at  the  contact  could  be 
approximately  applied  by  using  the  concept  of  "equivalent  radius".  This  equivalent 
radius,  Re  =  R1R2/  (Ri  +  R2),  where  Ri  and  R2  are  the  radii  of  the  particles  in 
contact,  has  been  used  in  the  past  in  the  case  of  unequal  cylinders  (Poritsky  1950). 
It  is  derived  from  the  Hertz  (1882)  theory  for  smooth  spheres  (f=0)  and  can  be 
applied  to  the  case  of  roug’i  spheres  as  a  first  approximation  (Ng  1989). 

The  above  array  was  first  subjected  to  isotropic  pressure,  a0  =  132  KPa;  the 
corresponding  geometrical  and  statistical  information  is  given  in  Figs.  39  and  40. 
As  the  spheres  now  are  not  identical  there  is  no  crystalization  and  the  531  sphere 
aggregate  is  almost  isotropic.  Therefore,  this  531  spheres  arrays  exhibits  a  higher 
level  of  symmetry  than  the  477  equal  spheres  arrays  and  can  be  assumed  to  be 
isotropic. 

The  array  was  then  loaded  under  biaxial  compression  to  a 22  —  332  KPa  while 
keeping  an  constant.  The  array  configuration  and  the  distribution  of  the  contact 
forces  at  the  end  of  loading  are  shown  in  Fig.  41;  the  corresponding  micromechanical 
information  is  shown  in  Fig.  42.  This  simulation  was  run  for  the  same  reasons  as 
with  the  477-sphere  array,  and  also  in  order  to  generalize  the  earlier  findings  for  the 
case  of  unequal  spheres. 

The  constrained  moduli  were  again  calculated  at  different  values  of  K  = 
auj <722-  The  normalized  constrained  moduli,  from  this  simulation,  as 

well  as  of  the  simulations  on  the  477  equal  sphere  array,  appear  in  Fig.  43,  where 
they  are  plotted  against  the  stress  ratio,  K  =  Figure  43a  contains  the 
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values  of  DAMP  spheres  in  the  direction  of  the  increasing  principal  stress  ( ^22), 
while  Fig.  43b  contains  the  normalized  moduli  in  the  other  direction,  in  which  the 
stress  (crn)  is  kept  constant.  Both  figures  also  include  data  points  from  a  number  of 
experiments  on  sand  (squares)  performed  at  the  University  of  Texas  by  Kopperman 
et  al  (1982),  as  well  as  the  results  of  the  nonlinear  finite  element  simulations  of 
Section  2.6. 

The  agreement  between  all  normalized  results  in  Fig.  43  is  quite  good  and  a 
straight  line  of  equation  =  (<^22/ <^ii)0' 38  =  K0'38  could  be  fitted  to  all 

numerical  and  experimental  data  points  included  in  Fig.  43a.  It  should  be  noted 
that,  if  the  change  in  sand  density  as  K  increases  is  neglected,  ~ 

[V£K)/Vp)]2.  Using  the  results  of  the  experiments  of  Kopperman  et  al  (1982)  in 
sand  one  computes  d£*P/d£P  ~  K0'4,  which  is  very  similar  to  the  line  in  Fig.  43. 
Therefore,  the  results  of  the  experiments  performed  at  the  University  of  Texas,  that 
the  P-wave  velocity  (Vp  =  jD/p)  propagating  along  a  principal  stress  direction  is 
in  first  order  approximation  only  a  function  of  that  principal  stress,  is  fully 
predicted  by  all  numerical  techniques  developed  as  part  of  this  research. 

To  further  demonstrate  the  fact  that  the  P— wave  propagation  velocity  (or  the 
constrained  modulus)  is  only  a  function  of  the  principal  stress  in  the  direction  of 
wave  propagation,  and  not  of  the  mean  stress  a0,  the  477  equal  sphere  array  was 
subjected  to  a  biaxial  compression-extension  loading  simulation  under  constant 
mean  stress.  a0  =  91  KPa.  The  vertical  stress  a 22 ,  was  increased  from  91  to 
137.4  KPa,  while  the  horizontal  stress,  an,  was  decreased  from  91  to  44.7  KPa. 
Thus  <70  =  0.5  (137.4  +  44.7)  =  91  KPa  at  the  end,  equal  to  the  initial  value  of  a0. 
Again  the  normalized  constrained  moduli,  M*P/d(P  were  computed  at  specific 
values  of  the  stress  ratio,  K  =  cth  /  a0  in  both  directions  an  and  an,  and  are 
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plotted  versus  the  stress  ratio,  K  in  Fig.  44.  If  the  moduli  were  a  function  of  the 
mean  stress,  the  P-wave  velocity  (and  modulus  in  that  direction)  would  be 
unaffected  by  the  changes  in  the  magnitude  of  the  principal  stresses  and  the  results 
would  plot  as  a  horizontal  straight  line  =  1.0.  This  is  not  the  case, 

however,  and  the  moduli  either  decrease  or  increase  together  with  the  change  in  the 
magnitude  of  stress  in  that  particular  direction.  Therefore,  Fig.  44  confirms  that 
the  constrained  moduli,  Dh,  are  affected  by  the  principal  stress  in  the  direction  of 
loading  only  and  are  not  by  the  mean  stress. 

Based  on  the  results  of  the  distinct  element  raicromechanical  numerical 
simulations  presented  herein,  it  is  concluded  that  changes  in  the  structure  of  the 
random  packings  and  granular  soils  are  responsible  for  the  observed  macroscopic 
behavior,  and  specifically  for  the  dependence  of  the  P— wave  velocity  only  on  the 
principal  stress  in  the  direction  of  wave  propagation.  Under  isotropic  compression 
(Figs.  36,  37,  39,  and  40),  the  orientations  of  the  contacts  have  a  more  or  less 
uniform  distribution,  the  material  behaves  isotropically  and  the  wave  propagation 
velocity  is  a  function  of  the  isotropic  stress,  cr0.  Under  biaxial  (2— D,  see  Figs.  41 
and  42)  or  triaxial  (3— D)  compression,  a  significant  number  of  contacts  is  gained  in 
the  direction  of  the  increased  major  principal  stress,  cri,  while  contacts  are  lost  in  all 
non— principal  directions.  These  contacts  along  principal  directions  form  "branches" 
or  "stiff  chains"  parallel  to  the  principal  stresses  which  transmit  most  of  the  applied 
principal  stresses  from  periodic  "boundary"  to  periodic  "boundary",  see  also  CundaJl 
and  Strack  (1983).  The  result  is  a  structure  reminiscent  of  a  simple  cubic  array 
loaded  by  stresses  <7\,  c-x  and  parallel  to  the  main  axes  of  the  array,  and  with  the 
effect  of  each  principal  stress  being  more  or  less  uncoupled  from  the  other  two.  This 
uncoupling,  which  in  the  simple  cubic  array  happens  naturally  due  to  the  geometry 
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of  the  array,  develops  in  the  random  arrays  as  a  consequence  of  these  stiff  columns 
or  "branches"  of  particles  which  carry  the  applied  load.  Therefore,  as  in  the  case  of 
the  simple  cubic  array  of  Section  2,  the  longitudinal  modulus,  D,  in  any  principal 
direction  is  a  function  only  of  the  stress  in  this  direction  and  is  unaffected  by 
variations  in  the  stresses  in  the  other  directions.  Although  the  CONBAL— 2 
simulations  discussed  herein  are  two  dimensional,  the  authors  postulate  that  the 
same  phenomenon  occurs  in  three— dimensions  with  the  formation  of  "stiff  chains"  in 
three  directions  instead  of  two.  It  has  been  shown  (Deresiewicz  1958,  1958a, 
Petrakis  and  Dobry  1986),  that  the  shear  moduli  of  the  simple  cubic  array  depend 
only  on  two  principal  stresses  and  are  independent  of  the  third  principal  stress. 
Since  the  moduli  of  the  random  arrays  are  expected  to  be  of  similar  nature  due  to 
the  3— D  "stiff  columns"  just  discussed,  it  is  believed  that  this  is  the  reason  why  the 
shear  moduli  (or  S— wave  velocities)  have  been  observed  experimentally  to  depend 
only  on  the  principal  stresses  in  the  direction  of  wave  propagation  and  particle 
motion.  However,  this  hypothesis  could  not  be  verified  with  the  2— D  simulation 
shown  and  further  work  needs  to  be  done  once  the  3— D  version  of  CONBAL  is 
developed. 

Despite  the  limitations  of  CONBAL— 2  discussed  above  (2— D  instead  of  3— D, 
no  particle  rotation),  the  results  presented  in  the  previous  section  and  by  Ng  (1989) 
show  that  the  program  is  a  powerful  tool  for  simulating  the  monotonic  and  cyclic 
behavior  of  sand.  Therefore,  it  will  also  be  used  to  investigate  the  validity  of  the 
assumptions  regarding  the  possible  distortion  of  the  shape  of  the  soil  yield  surface 
presented  in  Section  1.1. 


The  next  two  sections  present  a  combined  effort  including  laboratory  and 
numerical  experiments  aimed  at  defining  the  necessary  parameters  of  a  constitutive 
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law  for  granular  media,  while  Section  5  discusses  the  main  features  such  law  must 
have  in  order  to  reflect  the  observed  macroscopic  behavior.  These  parameters 
include  the  initial  and  subsequent  shape  of  the  yield  surface(s),  the  hardening  and 
flow  rules,  etc.  Both  the  laboratory  and  numerical  experiments  supplement  each 
other  as  they  correspond  to  similar  materials  and  stress  paths.  The  numerical 
simulations  provide  micromechanical  information  not  available  from  the 
experiments,  providing  additional  insight  and  helping  in  the  interpretation  of  the 
macroscopic  data. 


3.  EXPERIMENTAL  STUDY 


The  main  objective  of  this  experimental  work  is  to  determine  the  hardening 
characteristics  of  the  yield  surface(s)  of  a  granular  medium  after  experiencing 
various  load  histories.  For  this,  a  number  of  experiments  were  performed  along 
specific  stress  paths  on  hollow  cylindrical  specimens  composed  of  glass  beads.  These 
stress  paths  are  similar  to  those  used  by  Phillips  et  al.  (1972)  in  the  experimental 
study  of  the  yield  behavior  of  aluminum  discussed  in  Section  1.2.  Since  the  concept 
of  yield  surface  or  yield  function  is  central  to  plasticity  theory,  the  shapes  of  the 
initial  yield  surfaces  were  defined  and  measured  on  distinct  but  identical  specimens. 
As  a  second  step,  the  motion  and  the  shapes  of  these  yield  surfaces  after  loading 
were  also  obtained  by  a  series  of  cyclic  experiments  run  on  the  same  specimen. 

As  discussed  in  Section  1.3,  even  at  constant  pressure  soils  have  different  yield 
characteristics  than  metals,  with  the  most  important  difference  being  the  absence  of 
a  clearly  defined  elastic  region.  While  it  is  possible  to  define  an  initial  yield  surface 
through  the  threshold  strain  (for  example  the  locus  of  all  points  at  which 
irreversible  volumetric  or  pore  water  pressure  changes  start  to  occur),  this  is 
difficult  to  implement  in  the  laboratory.  An  alternative  is  to  define  the  yield 
surface  as  the  locus  of  all  points  having  a  certain  value  (close  to  the  threshold  but 
above  it)  of  the  octahedral  shear  strain;  this  approach  was  recently  adopted  by 
Peters  (1988)  and  will  also  be  used  in  this  study. 

The  experiments  were  performed  on  long  thin  hollow  cylindrical  specimens 
made  out  of  glass  beads.  Each  specimen  was  first  consolidated  isotropically  by 
having  the  same  internal  and  external  pressure,  cc{=  <r0),  and  then  was  subjected  to 
shear  bj  a  combination  of  axial  and  torsional  loads,  always  in  drained  condition, 
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using  a  computer— controlled  axial/torsional  servohydraulic  device.  In  any  given 
test,  the  mean  stress  was  kept  constant  and  equal  to  the  consolidation  pressure 
during  the  shear  stage.  Therefore,  all  experimental  data  from  the  test  are  on  one 
r- plane  in  stress  space. 


3.1  Laboratory  Equipment 

The  tests  were  performed  in  RPI’s  Class  of  1933  Earthquake  Engineering  and 
Cyclic  Loading  Laboratory.  This  laboratory  houses  the  MTS  axial/torsional  closed 
loop  servohydraulic  device  which  is  connected  to  a  digital  computer.  The  computer 
is  a  MACSYM-2  computer  made  by  Analog  Devices.  It  is  built  around  the  8086 
Intel  chip  and  uses  a  16-bit  bus  and  the  Intel  8087  co-processor  for  computing.  At 
the  front  end  it  has  a  8088  chip  for  fast  data  acquisition  and  control  from  ten 
available  channels.  Figure  45  sketches  the  configuration  of  this  computer  driven 
axial/torsional  apparatus. 

The  axial/ torsional  MTS  system  can  perform  cyclic  tests  in  both  axial  and 
torsional  modes,  alone  or  combined.  Each  mode  can  be  applied  either  monotonically 
or  cyclically,  in  stress  or  strain  control. 

The  applied  load  is  monitored  by  a  load  cell  which  can  detect  changes  in  the 
axial  force  of  ±  0.445X10 '2  KN.  The  axial  deformation  of  the  specimen  is  measured 
by  an  LVDT  located  at  the  top  of  the  loading  piston.  The  resolution  of  the  LVDT 
is  ±  0.635  X  10"4  m  Full  Scale  Output  (FSO). 
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The  applied  torque  is  measured  by  a  torsional  load  cell  which  is  able  to 
distinguish  changes  in  the  torque  in  the  order  of  ±  0.565  X  10 '4  KN— m.  The 
rotation  is  measured  by  an  RVDT  which  is  also  located  at  the  top  of  the  loading 
piston.  The  resolution  of  the  RVDT  is  ±  0.05  degrees  FSO. 

The  volumetric  changes  are  measured  by  an  WF17038  automatic  volume 
change  apparatus.  This  device  measures  the  changes  in  the  position  of  a  brass  piston 
incorporating  a  Bellofram  rolling  diagram.  Volumetric  changes  on  the  order  of 
±  0.05  cc  can  be  measured. 

The  cell  pressure  is  measured  by  a  pressure  transducer  which  is  connected  to 
the  cell  the  bottom  of  the  triaxi<u  :1.  The  resolution  of  the  pressure  transducer  is 
*  1.03  KPa  FSO. 

The  signals  from  the  previously  described  devices  are  transmitted  using  six 
different  channels  to  the  MACSYM-2  computer  through  a  Analog/ Digital 
converter;  its  resolution  for  the  particular  gain  used  in  these  tests,  was  *  4.88  mV. 
The  computer  performs  all  appropriate  calculations  and  sends  signals  to  keep  the 
mean  stress  constant  by  changing  the  cell  pressure.  This  is  achieved  by 
transmitting  the  signals  through  the  Digital/Analog  converter  to  Electropneumatic 
Regulator  NIT200.  This  regulator  is  capable  of  implementing  cell  pressure  changes 
of  ±  0.14  KPa  within  a  pressure  range  from  0  to  120  KPa. 


56 


3.2  Specimen  Properties 

Glass  beads  were  used  in  the  hollow  cylindrical  specimens.  Glass  beads  were 
used  instead  of  natural  sand  to  investigate  the  yield  surfaces  of  granular  media,  in 
order  to  minimize  the  variability  in  geometric  properties  which  may  influence  the 
results,  and  in  order  to  compare  with  the  numerical  simulations  using  program 
CONBAL— 2.  The  glass  beads  are  manufactured  by  Potters  Industries  Inc,  New 
Jersey,  and  are  made  from  soda— lime  glass. 

The  specific  gravity  of  the  glass  spheres  is  2.45  to  2.50  and  their  minimum 
roundness  is  70%.  The  Young’s  modulus  is  Es=68.97  X  106  KPa,  the  Shear 
Modulus  Gs=  29.66  X106  KPa,  and  the  Poisson’s  ratio,  us,  is  0.21.  The  static 
coefficient  of  friction  is  0.9  to  1.0  and  the  dynamic  coefficient  is  0.7  to  0.8.  The 
coefficient  of  friction  as  measured  by  Chen  (1986)  is  0.15.  Table  3  summarizes  the 
properties  of  the  material  used. 

The  hollow  cylindrical  specimens  had  an  outside  diameter,  D0,  of  7.11  cm  and 
an  inside  diameter,  Dj,  of  5.08  cm,  and  were  13.97  cm  tall,  as  shown  in  Fig.  46. 
These  dimensions  satisfy  the  criteria  for  stress  uniformity  in  such  specimens 
proposed  by  Wright  et  al.  (1978)  and  DeNatale  and  Vrymoed  (1981): 


Ri/R0  >  0.65 

(37) 

H  >  5.44  >/R3  -  R? 

(38) 
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where  Ri  and  R0  are  the  inside  and  outside  radii  of  the  hollow  cylinder,  respectively. 

All  specimens  were  prepared  by  the  dry  tamping  undercompaction  technique 
(Ladd  1978)  using  eight  layers  1.27  cm  thick.  The  undercompaction  value  of  the 
bottom  layer  was  chosen  equal  to  5%.  Undercompaction  was  used  instead  of  other 
alternative  techniques  because  it  produces  uniform  samples  which  can  be  easily 
reproduced  (Baziar  1987). 

Once  the  specimen  was  completed,  it  was  placed  in  the  triaxial  cell  which  was 
filled  with  deaired  water  to  prevent  migration  of  air  through  the  specimen 
membrane.  After  this  the  specimen  was  percolated  upwards  with  C02  to  eliminate 
air,  and  was  then  saturated  with  deaired  distilled  water  and  backpressured 
overnight  to  ensure  good  saturation. 

The  next  day  the  B  parameter  was  measured,  and  then  the  cell  with  the 
specimen  was  placed  in  the  MTS  axial/torsional  device. 

3.3  Testing  Program 

The  testing  program  consisted  of  a  series  of  monotonic  tests  and  of  another 
series  of  cyclic  tests.  The  monotonic  tests  were  used  to  classify  the  material  based 
on  its  monotonic  stress— strain  behavior  under  various  inclinations  of  principal 
stress,  while  the  cyclic  tests  were  used  to  define  the  yield  surfaces  (Kotsanopoulos 
1989). 

The  specimens  were  composed  of  a  mixture  of  sieve  40—50  (0.43  to  0.3  mm) 
and  sieve  60  —  80  (0.25  to  0.18  mm)  glass  beads  with  ratio  1:2  in  weight.  The  grain 
size  distribution  of  this  material  appears  in  Fig.  17.  A  mixture  was  used  instead  of 
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one  size  of  spheres  as  a  result  of  instabilities  and  sliding  which  occurred  during 
monotonic  tests  of  specimens  composed  of  identical  particles.  It  appears  that  these 
instabilities  and  sliding  are  the  result  of  regular  packings  and  crystalline  regions 
which  are  present  in  assemblages  of  equal  spheres,  but  not  necessarily  in  assemblies 
of  particles  of  unequal  diameters  All  specimens  were  prepared  by  the  same 
technique  described  in  Section  3.2  and  their  initial  void  ratio,  e,  was  the  same  in  all 
tests.  In  both  monotonic  and  cyclic  tests,  the  mean  stress,  <7C,  was  kept  constant 
during  loading  by  changing  the  cell  pressure  accordingly.  Therefore,  each 
experiment  corresponds  to  a  single  7r-plane  defined  by  this  ac. 

Five  drained  monotonic  radial  shear  tests  were  performed  as  listed  in  Table  3. 
In  three  of  these  tests,  the  specimens  were  isotropically  consolidated  to  ac  =  140 
KPa,  while  the  fourth  specimen  was  consolidated  to  277  Kpa  and  the  fifth  to  414 
KPa.  For  each  of  the  three  radial  shear  tests  performed  at  oc  £  140  KPa,  tne 
principal  stress  formed  a  constant  angle  u>  with  respect  to  the  vertical  direction, 
with  u  =  0°,  45°  and  90°.  These  correspond  to  vertical  compression,  torsional, 
shear,  and  vertical  extension,  respectively.  Ir  all  tests  the  mean  stress  was  kept 
constant  at  the  ac  of  the  test  (~  140  KPa).  The  other  two  tests  were  performed  in 
compression  (u  =0°).  The  test  conditions  and  stress  parameters  of  each  monotonic 
test  are  summarized  in  Table  4. 

The  two  last  columns  of  Table  4  include  the  angle  of  shear  direction,  0 ,  and 
the  coefficient  b  used  in  the  five  tests.  The  angle  0  defines  the  direction  of  the 
octahedral  shear  stress  in  the  octahedral  plane  and  is  given  by  tan  0  =  [3(<72  —  ff3)/ 
(2ai  -  <72  -  (73)|  '  .  The  coefficient  b  gives  the  relative  magnitude  of  the 
intermediate  principal  stress  a 2,  and  is  defined  as  b  =  (<72  -  o\)!(az  —  a\). 
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The  results  from  the  monotonic  tests  performed  at  ac  =  140  KPa  appear  in 
Figs.  48  and  49.  Figure  48  portrays  the  stress-strain  curves  (roct/flc  vs  70ct) 
obtained  from  these  tests;  it  can  be  seen  that  while  all  stress  paths  give  different 
peak  strengths,  the  shear  test  (ij=45°,  0=30°)  has  also  a  different  initial  stiffness; 
this  suggests  that  the  material  is  not  isotropic.  Figure  49  presents  the  measured 
volumetric  strain  as  a  function  of  the  octahedral  shear  strain,  and  here  again,  the 
behavior  is  different  depending  cn  the  stress  path  followed. 

The  resulting  stress— strain  curves  from  three  monotonic  compression  tests  at 
ac  =  136.8,  277.  and  414.1  KPa  appear  in  Fig.  50,  and  the  volumetric  behavior  is 
shown  in  Fig.  51.  Another  interesting  plot  is  the  variation  of  octahedral  shear  stress 
(roct)  with  volumetric  strain,  which  appears  in  Fig.  52;  plots  like  this  are  used  in 
constitutive  law  modelling,  since  they  indicate  at  which  stress  level  volumetric 
strain  starts  to  accumulate. 


3.4  Experimentally  Obtained  Yield  Surfaces 

As  discussed  previously,  granular  media  exhibit  nonlinear  inelastic 
stress-strain  behavior  even  at  very  small  strains.  Thus,  cohesionless  aggregates  do 
not  have  a  clearly  defined  elastic  region,  and  as  a  result  each  initial  yield  surface  can 
not  be  defined  the  way  it  is  done  in  metals;  as  the  locus  of  all  points  beyond  which 
plastic  deformation  occurs.  To  assume  a  mathematical  model  of  the  type  proposed 
by  Mroz  (1967)  or  Prevost  (1978),  which  describes  the  nonlinear  behavior  of 
materials  by  a  number  of  yield  surfaces  associated  with  a  constant  elastoplastic 
modulus  would  link  the  results  of  this  research  to  a  specific  type  of  model  and  would 
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also  be  very  difficult  to  achieve  in  the  laboratory.  Consequently,  a  yield  surface  was 
not  defined  here  as  the  locus  of  all  points  with  a  cmstant  elastoplastic  modulus. 
Instead,  the  yield  surface  was  defined  as  the  locus  of  all  points  in  stress  space  (rze, 
(<rzz  —  arr)/ 2  )  having  the  same  total  octahedral  shear  strain  70ct  =  0.05%  or  0.10% 
in  monotonic  tests.  This  definition  of  yield  surface  is  similar  to  that  used  by  Sture 
et  al.  (1987)  for  studying  the  yielding  characteristics  of  sand. 

Figure  53  sketches  the  initial  yield  surface  obtained  from  the  hollow 
cylindrical  tests  on  glass  beads  described  in  the  previous  section.  Each  point  of  the 
yield  surface  corresponds  to  a  different  drained  monotonic  radial  shear  test 
conducted  on  a  new  specimen. 

Then,  two  series  of  two  reversible  (cyclic)  experiments  each  were  conducted  in 
order  to  define  the  hardening  characteristics  of  the  yield  surface  when  the  specimen 
is  prestrained  in:  i)  compression,  and  ii)  shear  (torsion).  The  test  conditions  and 
stress  parameters  of  each  reversible  test  are  summarized  in  Table  5.  The  sequence 
of  loading  paths  for  these  tests  is  shown  in  Fig.  53,  where  the  numbers  may  be 
followed  to  identify  the  applied  stress  path.  Each  specimen  was  prestressed  to  point 
Rx  (x  =  a,  t)  in  stress  space  which  is  called  the  "Reference"  point.  Then  the  load 
was  increased  until  a  change  in  octahedral  strain  of  0.05%  was  detected  (points  1  in 
Fig.  53),  wherupon  the  load  was  reversed  to  Rx  and  then  further  decreased  to  point 
3,  defined  by  |  -  7^  |  =  0  05%.  This  probing  procedure  was  repeated  in 

different  directions  until  all  indicated  points  on  the  yield  surface  had  been  defined 
(see  Fig.  53).  Throughout  the  experiment,  the  mean  stress  was  kept  constant. 
Since  this  procedure  cannot  be  implemented  manually,  a  data  acqv.  ition/  control 
software  was  developed  to  control  the  above  experiments,  and  the  corresponding 
flowchart  is  shown  in  Fig.  54. 


61 


As  described  above,  stress  increments  such  as  used  above  (for  example  to  go 
from  points  1  to  3  in  Fig.  53)  were  either  a  torsional  or  a  compression  increment 
depending  on  the  test,  and  torsion  and  compression  increments  were  not  applied 
simultaneously,  so  as  not  to  distort  the  shape  of  the  yield  surface.  While  it  is  true 
that  in  order  to  reach  a  yield  surface  a  stress  must  be  applied  which  may  cause  the 
yield  surface  to  slightly  displace  or  distort,  the  stress  paths  described  above 
minimize  this  effect  by  cancelling  it  (e.g.,  stress  path  9—10—11—12  in  Fig.  53).  This 
occurs  because  the  yield  surface  is  expected  to  be  symmetric  about  either  the 
horizontal  deviatoric  stress  or  the  shear  axis,  depending  on  the  test.  This  procedure 
is  different  from  that  employed  by  Peters  (1988),  in  which  the  yield  surface  is 
determined  by  radial  stress  paths  on  different  specimens,  and  from  that  followed  by 
Phillips  et  al.  (1972),  see  Fig.  5.  On  the  other  hand,  it  is  very  similar  to  that 
utilized  for  composite  materials  by  Dvorak  et  al.  (1988). 

One  important  decision  in  this  work  was  to  use  the  prestress  point  Rx  as  the 
"reference"  point  from  which  all  stress  probes  start  Dvorak  et  al.  (1988)  used  the 
same  technique,  but  started  the  stress  probes  from  the  midpoint  between  points  1 
and  3  (Fig.  53)  in  stress  space.  Point  3,  as  well  as  all  other  points  on  the  yjeld 
surface  was  defined  by  Dvorak  et  al.  as  the  point  at  which  nonlinearity  starts.  Soils, 
however,  do  not  have  a  linear  region  beyond  which  permanent  deformations  occur, 
and  thus  another  criterion  must  be  used  to  select  the  reference  point,  which  is 
critical  because  it  controls  the  shape  of  the  yield  surface.  Prestress  point  Rx  was 
chosen  herein  for  reference  because  it  is  the  prestraining  to  this  point  which  controls 
the  subsequent  yielding  behavior  of  the  specimen.  Therefore,  it  was  decided  to  start 
all  stress  probes  from  point  Rx  in  stress  space.  If  soils  had  a  clearly  defined  linear 
elastic  and  plastic  regions,  this  method  would  coincide  with  that  of  Dvorak  et  al. 
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(1988).  Although  point  R*  remains  more  or  less  in  the  same  position  in  stress  space, 
there  is  permanent  deformation  associated  with  these  cyclic  stress  paths. 
Unfortunately,  this  is  a  result  of  the  properties  of  the  specimen,  and  little  can  be 
done  other  than  keeping  the  strains  as  small  as  possible.  On  the  other  hand,  if  a 
multiple  yield  surface  model  is  envisioned,  these  stress  paths  would  correspond  to 
the  stress  probes  needed  to  define  any  of  the  nested  yield  surfaces.  This  yield 
surface  would  then  be  analogous  to  those  of  Mroz  (1967)  or  Prevost  (1978),  since  it 
would  also  map  a  state  of  the  material  on  a  parameter,  in  this  case  70Ct- 

Figures  55  and  56  summarize  the  initial  and  subsequent  positions  of  yield 
surfaces  defined  from  monotonic  and  cyclic  tests  with  trc  ~  140  KPa.  A  total  of  four 
cyclic  tests  were  performed;  in  all  of  them  the  mean  stress  was  kept  constant  and 
equal  to  the  initial  confining  pressure,  ac  £  140  KPa.  Therefore  all  experimental 
data  lie  on  this  ir-plane.  Of  these  four  tests,  two  were  prestrained  in  shear  while 
two  were  prestrained  in  compression.  In  the  first  series  of  experiments,  the  yield 
surface  was  defined  as  the  locus  of  all  points  with  70Ct  =  0.1%.  In  ihe  second  series 
of  experiments  the  strain  level  defining  the  yield  surface  was  reduced  to  0.05%.  A 
further  reduction  of  the  strain  level  is  anticipated  in  the  future,  so  as  to  approach 
the  threshold  strain,  ^  1.0X10'4  to  2.0X10"4,  as  close  as  possible.  The  specimens 
of  the  second  series  of  the  experiments  were  also  subjected  to  600  cycles  of  torsional 
shear  strain  amplitude  7^  =  0.07  %  in  strain-control  under  drained  conditions  prior 
to  the  actual  experiment.  This  was  done  to  reduce  the  inherent  anisotropy  created 
during  the  sample  preparation  and  to  stiffen  the  sample  on  horizontal  planes. 
However,  it  was  found  that  this  neither  stiffened  nor  disturbed  the  specimen. 

The  results  of  the  first  series  of  tests  appear  in  Fig.  55.  The  yield  surfaces  are 
defined  as  the  loci  of  all  points  with  7oct=0.1%  and  were  obtained  by  prestraining 
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one  sample  in  torsion  to  7P  =  70Ct  =  0.3%  and  the  other  sample  in  compression  to 
7p  =  7oct  =  0.79%.  In  the  second  series  of  the  experiments  the  strain  yield  criterion 
was  reduced  by  50%,  to  7Oct~0.05%  (Fig.  56).  One  sample  was  prestrained  in  shear 
(torsion)  to  0.3%  and  a  second  sample  was  prestrained  in  compression  to  70ct  = 
0.5%.  In  both  series  of  tests  ac  was  kept  constant  at  about  1*0  KPa. 

The  results  suggest  a  significant  prestraining  effect  and  distortion  of  the 
chosen  yield  surface  in  the  direction  of  loading.  The  size  of  the  yield  surface  has 
decreased  in  size  in  the  direction  of  prestraining,  while  in  the  other  direction  the  size 
remains  essentially  unchanged.  Moreover,  it  is  observed  that  the  yield  surface 
translated  in  stress  space.  These  findings  are  important,  since  the  usual  hypothesis 
of  pure  kinematic  hardening  would  have  predicted  that  the  new,  translated  yield 
surface  size  would  be  the  same  as  before. 

Finally  a  "failure"  surface  was  constructed  out  of  all  monotonic  tests;  this 
surface  is  shown  in  Fig.  57,  where  the  deviatoric  stress  q=  is  plotted 

versus  p  =  +  cr3) /2.  The  shape  of  this  failure  surface  can  be  approximated  as 

conical.  It  is  reasonable  to  expect  that  this  failure  surface  must  be  very  similar  to 
the  yield  surfaces. 

These  two  series  of  cyclic  experiments  are  in  excellent  agreement  between 
themselves;  moreover,  they  are  in  general  qualitative  agreement  with  results 
obtained  in  other  materials  such  as  rr  jd  metal  matrix  composites.  Since  the 
common  plastic  deformation  mechanism  between  polycrystalline  aggregates  and 
granular  assemblies  is  sliding  taking  place  along  sliding  planes,  it  is  believed  that 
these  observed  results  are  a  consequence  of  those  sliding  mechanisms  occuring  at  the 
micromechanical  level. 
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4.  NUMERICAL  SIMULATIONS 

In  order  to  confirm  the  experimental  results  presented  in  Section  3  and  to  gain 
insight  into  the  micromechanical  processes  occuring  at  the  contact  level,  a  series  of 
2— D  numerical  simulations  was  performed  using  program  CONBAL— 2,  along  stress 
paths  very  similar  to  the  experimental  ones,  on  the  random  array  of  531  unequal 
spheres  presented  in  Section  2.7.  This  array  was  used  since  the  glass  beads  tested  in 
the  laboratory  was  composed  of  spheres  of  different  diameters. 

The  stress  paths  followed  are  very  similar  to  those  of  the  laboratory 
experiments  presented  in  Fig.  53,  with  two  minor  differences:  i)  the  octahedral 
strain  used  as  the  yield  criterion  was  first  set  to  7Oct=0.05%  and  was  later  reduced 
to  0.02%,  which  is  very  close  to  the  threshold  strain,  7t;  and  ii)  the  numerical 
simulations  were  strain  controlled,  thus  all  points  were  defined  in  strain  space. 
Because  of  the  two-dimensional  nature  of  the  numerical  simulations,  qualitative 
rather  than  quantitative  agreement  between  the  simulations  and  the  laboratory 
experiments  was  sought.  These  two  differences  between  numerical  simulations  and 
laboratory  experiment  are  not  very  important,  and  if  anything,  the  smaller  70Ct  used 
in  the  simulations  are  better,  since  a  smaller  70Ct  is  closer  to  the  threshold  strain 
and  thus  defines  a  region  near  the  cut-off  point  between  permanent  geometrical 
changes  and  effects  due  to  sliding  taking  place  at  the  contacts. 

The  531  quartz  array  was  consolidated  to  130  KPa;  the  particle  configuration 
and  the  micromechanical  statistics  at  this  stage  were  presented  in  Figs  39  and  40. 
Then,  two  series  of  calculations  simulating  four  tests  on  four  initial  identical 
specimen  were  conducted  in  which  the  array  was  prestrained  in:  i)  compression  to 
two  different  strain  levels,  and  ii)  in  shear  to  two  different  strain  levels.  All  four 
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"specimens"  were  prestrained  first  to  7oct=0.98%  and  subsequently  to  7OCt=2.33%. 
The  mean  stress  was  kept  constant  in  all  cases  at  ac=130  KPa  throughout  each 
simulation. 

Figure  5?  portrays  the  results  of  the  simulations  which  defined  the  yield 
surface  as  the  locus  of  all  points  with  7oct=0.05%.  There  are  five  yield  surfaces  in 
this  figure,  the  circular  one  being  the  initial  surface  which  was  defined  by  four 
monotonic  tests.  The  "subsequent"  yield  surfaces  correspond  to  prestraining  to  two 
strain  levels,  7P  =  70ct  =  0.98%  and  2.33%  as  marked.  It  can  be  seen  that  there  is  a 
significant  change  in  the  shape  of  the  prestrained  yield  surface:  it  has  shortened  in 
the  direction  of  loading  while  its  dimension  in  the  other  direction  has  not  changed. 
Moreover,  this  shrinking  of  the  yield  surface  increases  with  increasing  prestraining. 
Finally,  the  surface  became  flatter  in  the  direction  opposite  to  the  direction  of 
loading.  These  results  are  very  similar  to  those  obtained  in  the  laboratory 
experiments  on  glass  beads  (Figs.  55  and  56)  and  to  those  obtained  in  other 
materials  by  various  researchers  (Fig.  59). 

Figure  60  portrays  the  results  of  the  numerical  simulations  where  the 
magnitude  of  the  octahedral  shear  strain  defining  the  yield  surface  was 
70ct  =  0.02%.  The  results  are  again  very  similar  to  those  of  the  simulation  in  Fig. 
58,  to  the  laboratory  experiments  on  glass  beads  (Figs.  55  and  56)  and  on  other 
materials  (Fig.  59).  The  fact  that  the  yield  surfaces  defined  by  7Oct=0.02%  are 
quite  similar  in  size  to  those  defined  by  7Oct=0  05%  suggests  that  at  these  strain 
levels  the  granular  material  behaves  more  or  less  elastically.  Moreover,  it  has  been 
shown  (Ng  1989)  that  there  are  very  small  permanent  volumetric  changes  in  these 
numerical  simulations  at  a  strain  level  of  0.02%  or  less,  which  confirms  the  fact  that 
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the  effects  of  the  distortion  of  the  yield  surfaces  observed  is  not  a  consequence  of 
accumulated  strain,  but  rather  a  result  of  sliding  and  other  micromechanical 
phenomena  occurring  at  the  particle  level. 
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5.  PROPOSED  CONSTITUTIVE  LAW  FOR  GRANULAR  MEDIA 

This  section  presents  the  main  features  of  a  proposed  constitutive  law  for 
granular  media  using  the  results  presented  in  this  report.  Given  the  limitations  of 
the  existing  constitutive  laws  as  discussed  in  Section  1.1,  the  authors  have  taken  a 
different  approach  to  the  formulation  of  a  constitutive  law  for  dry  granular  soil. 
This  approach  has  as  a  starting  point  the  micromechanical  force— deformation 
contact  law  obtained  by  Seridi  and  Dobry  (1984)  and  Dobry  et  al.  (1989).  In  that 
contact  model,  there  is  an  infinite  number  of  yield  surfaces  of  conical  shape,  all 
parallel  to  each  other  at  all  times.  A  stress-strain  equivalent  of  this 
force-deformation  contact  law  will  be  used  for  the  proposed  constitutive  relation  of 
granular  media,  modified  to  incorporated  additional  important  features  observed  in 
granular  soils. 

The  model  is  based  on  the  results  of  the  laboratory  experiments  presented  in 
Section  3,  enhanced  by  the  numerical  simulations  discussed  in  Section  4.  The 
proposed  model  uses  the  most  important  features  of  the  contact  law  (Dobry  et  al. 
1989)  and  is  capable  of  incorporating  the  distortion  of  the  shape  of  the  yield  surfaces 
due  to  loading,  previously  discussed. 

Since  the  contact  model  was  naturally  formulated  in  terms  of  forces,  the 
mathematical  formulation  used  as  the  basis  of  the  proposed  constitutive  law  is  the 
viscoplasticity  model  presented  by  Yen  (1979)  and  Yen  and  Eisenberg  (1987),  but 
without  the  rate  effects.  This  model  must  be  modified  to  make  it  dependent  on  the 
mean  stress,  ac  =  p  =  <Ui/3,  and  the  hardening  and  flow  rules  must  be  changed  to 
match  those  of  the  contact  law  as  well  as  the  experiments.  The  yield  surfaces  are 
initially  conical,  and  they  are  permitted  to  distort  in  the  direction  of  loading  as 
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observed  in  the  experiments  and  numerical  simulations  (Figs.  55,  56,  58,  and  60). 
The  flow  rule  is  associative  only  on  the  deviatoric  plane  and  the  hardening  rule 
allows  for  parallel  translation  of  the  yield  surfaces  in  a  way  similar  to  that  of  the 
contact  model. 

The  main  aspects  of  the  proposed  constitutive  law  are  summarized  below.  In 
3— D  stress  space,  the  initial  deviatoric  stress  state,  Sij,  is  referred  to  the  subsequent 
state,  sy,  through: 

Sij  =  Sij  —  Qtij  +  Rij  (^9) 

where  Sij  is  the  deviator  stress,  ay  the  location  of  the  center  of  any  given  yield 
surface  (there  is  a  large  number  or  even  an  infinite  number  of  them)  in  deviatoric 
stress  space,  and  Ry  is  a  measure  of  the  distortion  of  this  yield  surface.  The  term 
Ry  does  not  appear  in  the  force-displacement  contact  model,  since  the  yield 
surfaces  were  not  allowed  to  distort;  it  is  used  here  for  the  first  time,  in  an  attempt 
to  present  a  model  which  is  general  enough  to  be  modified  to  take  into  account  the 
distortion  of  the  yield  surfaces.  Usually,  Ry  is  a  function  of  the  previous  deviatoric 
stress  state,  Sy,  and  of  the  history  of  material  deformation  (Eisenberg  and  Yen 
1984).  Assuming  a  yield  function  F  of  the  same  form  as  the  micromec^anical 
contact  law  (Dobry  et  al.  1989),  it  will  be: 

F  =  2  (sy  -  ay  +  R ij ) ( S i j  -  ay  +  Ry)  -  m2(p  -  pi)2  (40) 

where  m  is  a  material  constant  related  to  the  friction  angle  of  the  soil  and  defined 
through  laboratory  experiments,  p  =  ^  0-51  is  the  value  of  the  hydrostatic  stress  for  a 
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point  on  the  cone,  while  pi  is  the  hydrostatic  stress  at  the  apex  of  the  same  cone. 
Figure  61  presents  several  conical  yield  surfaces  and  a  deviatoric  plane  (x  —  plane) 
in  principal  stress  space.  Figure  62a  shows  the  intersection  of  one  of  these  yield 
functions  with  ax  —  plane,  which  portrays  the  simultaneous  translation  (kinematic 
hardening)  and  distortion  of  the  resulting  yield  circle.  Figure  62b  defines  the 
corresponding  distortion  term  Ry  used  in  Eq.  40. 

If  the  distortion  term,  Rij  is  eliminated,  the  yield  function  reduces  to 

F  =  \  (sij  -  <*ij)(sij  -  orij)  -  m2  (p  -  Pi)2  (41) 

which  is  the  exact  stress  equivalent  of  the  micromechanical  force— displacement 
contact  law  (Dobry  et  al.  1989). 

The  flow  rule  is  described  by: 

=  2~ffp  nijKnkl  +  umni)  (42) 

where  Hp  is  the  plastic  tangent  modulus.  In  Eq.  42,  njj  is  the  deviatoric  part  of  the 

!! 

normal  vector  to  the  yield  surface,  njj,  and  u^m  is  the  spherical  part  of  another 
vector  Uki,  so  that  normality  is  observed  only  on  the  deviatoric  plane.  An  example 

tf 

of  the  form  of  Umm  is  given  by  Prevost  (1985): 

=  (43) 

6  ivh)  +  i 
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3  1/2 

where  77  =  ^  (sy  Sy)  /p  is  a  stress  ratio,  and  rj  is  a  material  parameter.  Therefore, 
this  observation  of  normality  only  on  the  deviatoric  plane  is  a  characteristic  of  the 
micromechanical  constitutive  law,  as  well  as  of  a  number  of  other  models. 

The  hardening  and  distortion  parameters,  ay  and  Ry,  can  be  either 
determined  a  priori  or  from  a  combination  of  results  of  micromechanical  simulations 
and  laboratory  experiments.  The  exact  form  of  these  parameters  is  not  currently 
known  for  the  case  of  sand,  but  without  loss  of  generality  it  can  be  assumed  that 


ay  —  /ivy 

(44) 

Ry  =  — Auy 

(45) 

where  vjj  and  Uij  are  tensors  defining  the  specific  direction  of  hardening,  and  (i  and 
A  are  scalars  defining  the  magnitude  of  the  hardening  mechanism.  If  a  finite 
number  of  yield  surfaces  is  used  in  the  model,  the  kinematic  hardening  direction, 
vy,  can  be  defined  in  a  number  of  ways  including  those  proposed  oy  Prager  (1955), 
Ziegler  (1958),  Mroz  (1967),  and  Phillips  and  Weng  (1975).  However,  if  functions 
defining  an  infinite  number  of  yield  surfaces  are  used  such  as  done  in  the  contact 
law,  no  such  rule  is  necessary.  The  two  scalars,  n  and  A,  are  related  to  each  other 
through  the  consistency  condition. 

Yen  (1979)  and  Yen  and  Eisenberg  (1987),  using  the  results  from  experiments 
by  Phillips  and  his  coworkers  have  defined  both  ay  and  Ry  and  used  them 
successfully  in  a  constitutive  law  for  aluminum.  The  excellent  agreement  between 
theory  and  experiment  in  the  case  of  cyclic  loading  of  aluminum  is  shown  in  Fig.  63, 
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where  the  discrete  points  correspond  to  the  laboratory  data  and  the  solid  lines  to  the 
model. 

As  already  noted,  the  above  proposed  stress  equivalent  model  to  the  contact 
law  summarized  by  Eqs.  39-45  has  a  number  of  advantages  over  the  existing 
plasticity  models  for  sand,  while  still  lacking  in  some  respects.  The  advantages 
include  its  foundation  on  basic  micromechanical  features  such  as  the  behavior  at  the 
interparticle  contacts,  and  its  ability  to  model  a  number  of  important  observed 
macroscopic  features  such  as  the  pressure  dependence  of  stiffness,  the  generalized 
Mohr— Coulomb  failure  law  and  the  basic  hardening  behavior  during  cyclic  shear 
loading.  Also,  as  already  mentioned  and  included  in  Eq.  40,  the  model  can 
incorporate  the  distortion  of  yield  surfaces  due  to  prestraining.  Moreover,  the 
micromechanical  basis  of  the  model  allows  it  to  be  compared  directly  with  the 
contact  model  through  numerical  simulations  on  random  packings  of  spheres,  so 
that  the  model  can  be  adjusted  at  the  macroscopic  level  with  due  consideration  to 
the  underlying  micromechanical  phenomena. 

On  the  other  hand,  in  its  present  state  the  basic  model  of  Eqs.  39-45  still 
misses  some  features  which  have  been  observed  in  sands.  The  proposed  law  does  not 
model  the  dilation  ',T  ;ved  in  medium  and  dense  granular  soils  for  shear  loading 
contained  in  the  ir— plane,  and  it  does  not  account  for  inelastic  volumetric  strains 
due  to  isotropic  consolidation,  nor  for  overconsolidation  effects  Finally,  it  does  not 
consider  initial  (inherent,  structural,  elastic)  anisotropy,  which  would  permit  the 
prediction  of  the  anisotropic  wave  propagation  phenomena  reported  by  Stokoe  and 
his  co-workers  for  anisotropically  consolidated  sand  (Knox  et  al.  1982,  Kopperman 
et  al.  1982). 
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The  dilation  can  be  handled,  in  principle, by  modifying  both  the  hardening  and 
flow  rules.  Existing  models  such  as  Prevost’s  (1985)  account  for  this  by  allowing 
rotation  as  well  as  translation  of  the  yield  surfaces  (Fig.  64).  To  account  for  this, 
the  yield  surfaces  in  our  proposed  model  could  be  permitted  to  rotate  after  a  certain 
strain  level  has  been  reached,  whereupon  dilation  will  occur  (Fig.  65).  This  is 
consistent  with  the  experimental  finding  in  sand  that  dilation  occurs  in  medium  and 
dense  sand  only  after  a  certain  strain  level  beyond  the  threshold  has  been  reached. 
On  the  other  hand,  the  experimental  evidence  presented  in  Fig.  52  as  well  as  other 
observations  on  sand  from  the  literature  suggest  that  dilation  under  shear  occurs  at 
relatively  large  strains.  Therefore,  the  flow  rule  does  not  have  to  be  modified  at 
small  strains. 

The  inelastic  volumetric  strains  during  isotropic  consolidation  and  the 
response  of  overconsoli dated  soil  can  be  obtained  with  the  help  of  a  "cap"  (Roscoe  et 
al.  1958,  DiMaggio  and  Sandler  1971,  Baladi  and  Rohani  1979). 

The  problem  of  inherent  (structural)  anisotropy  —  especially  important  for 
anisotropically  consolidated  sand  —  is  quite  complex  even  in  the  case  of  metal 
plasticity.  "Inherent  anisotropy"  means  that  the  material  is  anisotropic  in  its 
reference  state.  After  plastic  flow  has  occurred,  in  general,  the  material  ceases  to 
have  the  anisotropy  of  the  reference  state.  A  very  good  example  is  the  case  of  rolled 
steel;  in  its  stress  free  state,  before  rolling,  the  material  is  isotropic  and  its  behavior 
could  be  described  by  isotropic  functions.  After  the  sheet  of  steel  has  been  subjected 
to  rolling,  it  may  exhibit  orthotropic  symmetry,  and  if  this  is  to  be  the  reference 
state,  the  material  behavior  should  be  described  by  orthotropic  functions  with 
respect  to  the  new  reference  state.  The  case  of  soil  anisotropy  is  very  similar:  if  the 
behavior  of  an  anisotropically  consolidated  sand  is  to  be  described,  this  has  to  be 
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done  with  orthotropic  functions.  In  general,  simple  solutions  applied  only  to  the 
elastic  part  of  the  strain,  such  as  proposed  by  Hardin  (1980),  are  not  appropriate. 
For  example,  if  the  tests  of  Stokoe  and  his  coworkers  are  to  be  taken  into  account, 
orthotropic  functions  are  needed  for  the  description  of  the  phenomenon.  It  is 
reported  by  Hill  (1950)  and  Dafalias  (1979)  that  the  yield  condition  must  take  into 
account  the  material  symmetries  present  in  the  initial  (reference)  state.  The 
direction  of  flow  is  dependent  upon  those,  and  the  yield  surface  must  take  them  into 
account.  This  not  only  complicates  the  expressions  of  the  yield  surface 
considerably,  but  the  implementation  of  the  law  as  well.  For  example,  according  to 
Dafalias  (1979),  the  expression  for  the  yield  condition  of  an  orthotropic  material 
becomes: 

f  =  HjjkiSijSki  —  k2  (46) 

where  Hijki  is  a  fourth  rank  tensor  which  depends  on  the  material  symmetry  and 
the  plastic  strain  and  k  is  the  yield  stress.  Despite  its  obvious  complexity,  the 
plasticity  of  initially  anisotropic  materials  can  be  taken  into  account  in  a  number  of 
ways  for  the  cases  of  transverse  isotropy  and  orthotropy.  The  case  of  transverse 
isotropy  is  very  common  in  the  mechanics  of  composite  materials  (Dvorak  and 
Bahei— El  Din,  1982),  and  is  handled  with  the  aid  of  four  stress  "pseudo— invariants" 
which  take  into  account  the  symmetries  of  the  material.  In  the  case  of  soils, 
Dafalias  (1986),  Dafalias  and  Herrmann  (1986)  and  Anandarajah  and  Dafalias 
(1986)  have  proposed  a  general  plasticity  model  which  has  been  applied  to  model 
isotropic  and  anisotropic  clay  in  undrained  condition.  The  case  of  drained  loading 
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of  sand  is  more  complex,  since  its  behavior  is  pressure  dependent  and  its  material 
symmetry  is  affected  by  variations  in  the  mean  stress  uii/3. 
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6.  CONCLUSIONS 

The  results  of  a  systematic  micromechanical  study  was  presented  on  the 
stress-strain  behavior  of  granular  media,  aimed  at  developing  a  constitutive  law  for 
sands  and  other  granular  materials. 

As  summarized  in  Table  1,  the  research  focused  on  the  computation  of  the 
stress  strain  response  of  arrays  of  rough,  elastic  quartz  spheres  of  similar  sizes 
simulating  uniform,  rounded  quartz  sand.  It  started  from  the  RPI  general  solution 
for  the  nonlinear  force  displacement  relation  at  the  contact  between  two  spheres 
(program  CONTACT),  and  it  continued  with  analyses  of  regular  and  random  arrays 
of  spheres  using  a  variety  of  analytical  and  numerical  techniques.  The  research 
progressed  from  very  small  (7  2  KH)  to  very  large  strains  (7  a  10'2  to  10'1),  and  it 
included  comp  risons  at  different  strain  levels  between  the  results  obtained  with 
these  techniques. 

A  number  of  regular  arrays  of  spheres  were  analyzed,  and  then  they  were 
combined  into  regular /random  arrays  to  attain  arbitrary  macroscopic  void  ratios. 
These  regular /random  array  models  were  loaded  isotropically  and  their  macroscopic 
response  was  evaluated  at  very  small  strains  by  the  Self  Consistent  Method.  Then, 
another  regular/ random  array  model  was  developed  and  analyzed  by  a  2D  Nonlinear 
Finite  Element  method.  This  Finite  Element  model  was  loaded  first  isotropically 
and  then  anisotropically  to  "failure"  (corresponding  to  the  beginning  of  generalized 
sliding  and  particle  rearrangement),  and  it  provided  information  from  7  «  10 ‘fl  up  to 
7  s  10 '4.  Finally,  2D  numerical  simulations  of  random  arrays  were  implemented 
using  a  Nonlinea-  Distinct  Element  technique,  covering  the  whole  strain  range  from 
7  »  10 to  10 _1.  In  addition,  comparisons  between  these  calculations  and 
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measurements  in  sands  and  other  granular  materials  were  performed,  including 
experiments  available  from  the  literature  and  axial-torsional  tests  on  glass  beads 
conducted  as  part  of  this  research. 

It  was  found  that  all  computational  techniques  are  in  very  good  agreement 
between  themselves  as  well  as  with  the  experimental  results.  Some  specific 
conclusions  for  the  different  strain  levels  are  as  follows: 

1)  In  a  granular  medium  isotropically  loaded  under  a  given  macroscopic 
pressure,  there  is  a  unique  relation  between  the  shear  stiffness  at  very  small  strains 
(Gmax)  and  the  average  number  of  load-transmitting  contacts  per  particle,  with 
Graax  increasing  with  the  number  of  contacts.  This  was  confirmed  by  all  techniques 
used  in  the  report. 

2)  In  the  same  isotropically  loaded  granular  medium,  the  relation  between 
void  ratio  and  Gm3r  is  not  unique.  For  the  same  void  ratio  and  particle  geometric 
arrangement,  Gmax  can  vary  by  a  factor  of  three  depending  on  the  number  of 
contacts.  In  a  typical  sand,  Gmax  is  much  smaller  than  predicted  from  a  regular 
array  or  combination  of  regular  arrays  having  the  same  macroscopic  void  ratio  as 
the  sand.  However,  millions  of  cycles  of  small  strain  amplitude  (10*4  <  7  <  10 _3) 
will  increase  the  stiffness  of  the  sand  and  make  it  approach  the  theoretical  value 
predicted  by  the  Self  Consistent  method,  by  increasing  the  number  of  contacts  while 
the  void  ratio  stays  essentially  constant. 

3)  The  P-wave  velocity  through  an  anisotropically  loaded  2D  granular 
medium  depends  mainly  or  exclusively  on  the  principal  stress  in  the  direction  of 
wave  propagation.  This  seems  to  be  the  result  of  a  specific  force  transmission 
mechanism  throughout  the  medium,  with  the  load  paths  oriented  more  or  less 
parallel  to  the  directions  of  the  applied  principal  stresses.  In  this  way,  the  behavior 
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of  the  medium  resembles  that  of  a  simple  cubic  array,  where  the  stiffness  in  a 
principal  direction  is  not  influenced  by  the  stresses  in  the  other  direction.  This 
result  was  confirmed  for  2D  by  both  the  Nonlinear  Finite  Element  model  and  the 
Nonlinear  Distinct  Element  simulations.  Based  on  this  analysis  and  on  available 
experimental  results  on  sand,  it  is  hypothesized  that  a  similar  conclusion  and  force 
transmission  mechanism  is  also  valid  in  3D,  for  both  P-wave  and  S-wave  velocities. 

4)  In  the  small  strain  range  below  the  threshold  (7  <  7t  s  10'4  for  quartz 
spheres),  the  shear  stress— strain  behavior  of  the  medium  is  controlled  by  the 
original  arrangement  of  the  grains  and  by  the  nonlinear  force— deformation  response 
at  the  contacts,  and  essentially  no  rearrangement  of  particles  takes  place.  Aspects 
of  this  small  strain  behavior  including  the  existence  of  the  threshold  were  known 
from  previous  experimental  and  analytical  research.  The  Nonlinear  Finite  Element 
model  confirmed  this  behavior.  It  also  provided  additional  insight  on  the  mechanics 
of  small  strain  response,  and  on  the  process  which  culminates  in  particle  separation 
and  rearrangement  when  the  shear  strain  reaches  the  threshold. 

5)  A  major  role  in  the  shear  stress-strain  behavior  of  granular  media  is 
played  by  sliding  taking  place  between  the  grains.  Available  analytical  and 
experimental  studies  in  polycrystalline  aggregates  such  as  metals  suggest  that 
specific  sliding  directions  at  the  interparticle  contacts  control  the  features  of  specific 
portions  of  the  yield  surface  of  the  medium  in  stress  ^  ce.  This  was  verified  with 
our  Nonlinear  Finite  Element  calculations  of  random/regular  arrays  of  spheres.  One 
corollary  of  these  studies  in  metals  is  that  the  yield  surface  not  only  translates  due 
to  shear  loading  but  also  shrinks  and  distorts  in  the  direction  of  loading,  forming  a 
smooth  vertex  at  the  loading  point.  This  same  behavior  was  observed  in  our  2D 
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simulations  of  random  arrays  of  quartz  spheres  using  the  Nonlinear  Distinct 
Element,  as  well  as  in  our  axial-torsional  tests  on  glass  beads. 

6)  A  constitutive  law  for  granular  media  was  proposed  based  on  these 
findings,  which  is  the  stress— strain  equivalent  of  the  RPI  force  deformation  model 
for  the  interparticle  contact,  enhanced  to  incorporate  dilation  for  strains  above  the 
threshold,  and  to  allow  for  distortion  of  yield  surfaces  in  the  direction  of  loading. 
Basic  characteristics  of  this  model  include  an  infinite  number  of  originally  parallel 
yield  surfaces  of  conical  shape,  as  well  as  modified  normality  and  kinematic  strain 
hardening  rules. 
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TABLE  2.  RPI  COMPUTER  PROGRAM  CONBAL-2 

(after  Ng  and  Dobry  1989) 

CQNBAL-2  =  CQNTact  truBAL  in  2-D 


TRUBAL  (Peter  Cundall) 

CONTACT  (RPI) 

Distinct  Element  Method 

General  Solution  to  Mindlin-Hertz 

Contact  Problem; 

3D  Random  Array  of  Spheres 

Reproduces  identically  all  Available 
Analytical  Results; 

Periodic  Space 

Program  Accepts  as  Input  Arbitrary 

Time  History  of  Force  or  Displacement; 

Linear  Pressure-Independent 

Can  be  Used  Approximately  for 

Arbitrary  Contact  Compliance 

Spheres  of  Somewhat  Different  Size 

I 


I 


CONBAXi-2 


TRUBAL  +  CONTACT 


ID  +  VECTORI ZATION 


I 

I 

I 

I 

I 


TABLE  3.  Properties  of  Glass  Beads 


Specific  Gravity 

Roundness 

Poisson’s  Ratio 

Young’s  Modulus 

Rigidity  (Shear)  Modulus 

Static  Coefficient  of  Friction 

Dynamic  Coefficient  of  Friction 


2.45  to  2.50 
Minimum  of  70% 
0.21 

68.97  x  10«  KPa 
29.66  x  108  KPa 
0.9  to  1.0 


0.7  to  0.8 


TABLE  4.  Summary  of  Monotonic  Radial  Shear  Tests  on  Glass  Beads  Specimens 


Type 

of 

Test 

(w°) 

Isotropic 

Confining 

Pressure 

(KPa) 

Void 

Ratio 

e 

Relative 

Density 

Dr 

% 

Pore 

Pres. 

Param. 

B 

Angle 

of 

Shear 

Direct. 

& 0 

Coefficient 

b 

Compr.  (0°) 

136.8 

0.593 

53 

0.954 

0 

0.0 

Tors. (45°) 

141.2 

0.603 

47 

0.963 

30 

0.5 

Ext.  (90°) 

140.0 

0.602 

48 

0.985 

60 

1.0 

Compr. (0°) 

277.0 

0.605 

46 

0.965 

0 

0.0 

Compr.  (0°) 

414.1 

0.603 

47 

0.955 

0 

0.0 

TABLE  5.  Summary  of  Reversible  Tests  on  Glass  Bead  Specimens 


Direction 
of  Pre- 
Straining 

Isotropic 

Confining 

Pressure 

(KPa) 

Void 

Ratio 

e 

Relative 

Density 

Dr 

% 

Pore 

Press. 

Param. 

B 

P  re¬ 
cycling 

Yield 

Surface 

Assoc. 

with 

1 A  7oct  | 

Compres. 

137.0 

0.595 

52 

0.987 

Yes 

0.05% 

Torsion 

138.1 

0.598 

50 

0.950 

Yes 

0.05% 

Compres. 

139.2 

0.597 

51 

0.954 

No 

0.1% 

Torsion 

139.5 

0.604 

47 

0.978 

No 

0.2% 

(*)  Applied  cyclic  shear  strain  7^  =  0.07%.  Number  of  cycles  N  =  600 
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Effect  of  stress  inauced  anisotropy  on  compressions!  wave  velocity  for 
triaxial  confinement  (Kopperman  et  al.  1982). 


Shear  Modulus,  G  max  (ksi) 


Cycles  of  High  Amplitude  Prestrain 


Figure  2.  Small  Strain  Shear  Modulus,  GBax,  versus  Cycles  of  Shear  Prestrain  — 
Hollow  Cylindrical  Samples  (Drnevich  1967). 


Cubical  Triaxial  Cell 
Experiments 

Oi 


Figure  3. 


Experimentally  obtained  yield  surface  of  sand  from  cubical  triaxial 
experiments  (Peters  1988). 
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Loading  paths  for  (a)  initial, 
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Experimentally  obtained  initial  and  subsequent  yield  surfaces  for  aluminum 
at  70°F  temperature  (Phillips  and  Tang  1972). 


Simple  Cubic  Unit  Cell  Simple  cubic  structure 


Figure  7. 


Simple  cubic  array  of  spheres. 
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Figure  8a. 


Elastic  Spheres  Under  Normal  and  Tangential  Loads. 
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Figure  8b.  Tangential  Force— Displacement  Relation  for  N  Constant, 
(Dobry  et  al.  1982). 
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Figure  10.  Force-deformation  Behavior  of  Two  Elastic  Rough  Spheres  in  Contact: 

Analytical  Solution  for  Oscillating  Oblique  Forces  with  (a)  dT/dN  >  f  and 
(b)  dT/dN  <  f;  (c)  N  Increasing,  T  Increasing;  (d)  N  Decreasing,  T 
Increasing;  (e)  N  Increasing,  T  Decreasing;  (f)  N  Decreasing,  T  Decreasing; 
(Mindlin  and  Deresiewicz  1953). 


TANGENTIAL  FORCE  (KG) 


Figure  11. 


TANGENTIAL  DISPLACEMENT,  <5(CM)(x10'4) 


Numerical  Simulation  of  the  Force— Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load-Displacement  Relation  for  an 
Oscillating  Oblique  Force  with  dT/dN  >  f. 


TANGENTIAL  FORCE  (KG) 


Figure  12. 


Numerical  Simulation  of  the  Force-Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load-Displacement  Relation  for  an 
Oscillating  Oblique  Force  with  dT/dN  <  f. 
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Figure  13. 


Numerical  Simulation  of  the  Force-Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load— Displacement  Relation  for  the 
Case  in  which  N  Increases  and  T  Increases. 


TANGENTIAL  DISPLACEMENT,  <5(CM)(x10‘3) 


Numerical  Simulation  of  the  Force-Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load— Displacement  Relation  for  the 
Case  in  which  N  Decreases  and  T  Increases. 
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Figure  15.  Numerical  Simulation  of  the  Force— Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load— Displacement  Relation  for  the 
Case  in  which  N  Increases  and  T  Decreases. 
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TANGENTIAL  DISPLACEMENT,  <5(CM)(x10'3) 

Figure  16. 


Numerical  Simulation  of  the  Force— Deformation  Behavior  for  Two  Equal, 
Elastic,  Rough  Spheres  in  Contact:  Load— Displacement  Relation  for  the 
Case  in  which  N  Decreases  and  T  Decreases. 
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Simple  cubic  structure 


Body-centered  cubic  unit  cell  Body-centered  cubic  structure 


Face-centered  cubic  unit  cell  Face-centered  cubic  structure 


Figure  17.  Three  Regular  Arrays  of  Equal  Spheres. 
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Figure  19. 


Shear  Wave  Velocity  versus  Void  Ratio:  Comparison  between  Predictions 
from  Regular  Cubic  Arrays  of  Quartz  Spheres  and  Measurements  in  Quartz 
Sand. 


Velocity  (ft/sec) 


o  High  Tolerance  Spheres  ( <p  =  V»  ±  10  *  10'6) 
*  Low  Tolerance  Spheres  (<p  =  V«  ±  50  *  I0~fl) 


Figure  20. 


Rod  Wave  Velocity  Measurements  in  Regular  Dense  Array  of  Steel  Spheres 
Loaded  Isotropically  (Duffy  and  Mindlin  1957). 


Confining  Pressure,  ao  (psi) 
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I 

I  Figure  23.  Isotropic  Confining  Pressure  versus  Volumetric  Strain:  Analytical  and 

®  Experimental  Results  for  e  =  0.54. 
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Figure  24. 
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(a)  Small  Strain  Shear  Modulus,  GBax,  and  (b)  Void  Ratio,  versus  1 Number 
of  Cycles  of  Prestrain  -  Solid  Cylindrical  Samples  (Song  and  Stokoe  1987). 
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Figure  25.  Shear  Modulus,  GmiX,  versus  Coordination  Number  for:  (a)  Regular  Arrays 
and  (b)  Random  Arrays  of  a  Given  Porosity  by  the  Self-Consistent  Method 
and  the  Analytical  Expressions  of  Walton  (1987).  All  Spheres  have  been 
Assigned  the  Properties  of  Quartz. 
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Figure  26. 


Stress  State  on  a  Simple  Cubic  Array  of  Identical  Spheres. 


TANGENTIAL  DISPLACEMENT  IN  cm  (x  10'2) 


TANGENTIAL  DISPLACEMENT  IN  cm  (x  10‘2) 


SHEAR  STRAIN  (x  10'2) 


Stress— Strain  Behavior,  (c),  of  a  simple  Cubic  Array  of  Quartz  Spheres 
Subjected  to  Biaxial  Stress,  <7U  =  161,  =  200  KPa,  Followed  by  Pure 

Shear,  <ru,  to  failure.  (a)  and  (b)  Portray  the  Corresponding 
Force— Deformation  Behavior  of  the  "weak1'  and  "strong"  contact, 
respectively. 


Figure  27. 


Figure  29. 


(a)  Orientation  of  the  Applied  Principal  Stresses  During  Pure  Shear 

Loading  with  Mean  Stress  Constant,  erg  =  0.5(cr^  +  cr$)  =  constant,  and 

(b)  orientation  of  an  element. 
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Figure  30.  Difference  in  Applied  Principal  Stresses,  crt  -  o-j,  versus  Shear  Strain, 
7  =  €i  —  cj,  Calculated  for  all  Media  (16  and  64  elements),  and  for  Values 
of  a  =  0°,  22.5°,  45°,  67.5°,  and  90°  (MX  =  medium  X,  A YT  =  inclination 
of  <T\  is  YY  degrees). 


HYDROSTATIC  PRESSURE,  al  (KPa) 
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Figure  31. 


Medium  1:  Stress— Strain  Behavior  for  Hydrostatic  Compression  Loading 
Starting  from  a*  =  300  KPa. 
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Figure  32.  Medium  1:  Stress-Strain  Behavior  for  Cyclic  Isotropic  Loading,  ag  ±  Aog, 
with  og  =  300  KPa  and  A<jg  =  200  KPa. 
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SHEAR  STRAIN  (x  10‘3) 


Medium  1:  Stress— Strain  Behavior  for  Cyclic  Pure  Shear  Loading  Following 
Isotropic  Compression,  <r<>  =  100  KPa.  (a)  7C  =  40  KPa  and,  (b)  7c  =  43 

KPa. 
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Figure  34. 


SHEAR  STRAIN,  y,  PERCENT 


Normalized  Shear  Modulus,  G/Gaax,  (a),  and  Damping  ratio  £,  (b),  versus 
Shear  Strain,  7.  Comparison  Between  Values  Measured  in  Sand,  the 
Simulated  Aggregate  and  Simple  Cubic  Array  in  Pure  Shear. 


Normalized  Constrained  Modulus, 


Figure  35. 


Normalized  Constrained  Moduli,  dPPMP  versus  Stress  Ratio, 
K  =  ay  Comparison  Between  Analytical  and  Experimental  Results  (a) 
in  the  Direction  of  o*l,  and  (b)  in  the  Direction  of  a%  is  kept  constant). 


°i  2  =  91  KPa 


Figure  36.  2-D  Random  Array  of  477  Equal,  Elastic,  Rough,  Quartz  Spheres 

Subjected  to  Isotropic  Compression,  <7  =  91  KPa.  Note  that  this  Figure 
Represents  the  "Window"  with  the  Representative  Random  Pattern. 
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Figure  37.  Statistical  Information  Regarding  the  isotropic  Compression  at  a0  -  91 
KPa  of  the  477-sphere  medium  of  Figure  36. 
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Figure  38. 


Shear  Modulus,  GBax)  vs  Coordination  Number  for  :  i)  Two  Random  Arrays 
of  477  Equal  Spheres,  ii)  Regular  Arrays,  and  iii)  Random  Arrays  of  a  given 
porosity  by  the  Self  Consistent  Method  and  the  analytical  expressions  of 
Walton  (1987). 
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Figure  39. 


2-D  Random  Array  of  531  Elastic,  Quartz  Spheres  of  Two  Radii 
(Ri/Ra  =  1.5)  Subjected  to  Isotropic  Compression  <r0  =  132  KPa. 
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Figure  40. 


Statistical  Information  for  the  isotropic  loading  of  Fig.  39. 
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Figure  41. 


2-D  Random  Array  of  531  Elastic,  Quartz  Spheres  of  Two  Radii 
(R1/R2  =  1.5)  Subjected  to  Biaxial  Compression  Loading. 
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Figure  42. 


Statistical  information  for  the  biaxial  compression  loading  in  Fig.  41. 


Figure  43. 


Normalized  constrained  moduli,  Dy i '  /  D\i',  vs  stress  ratio,  K  =  crj  /  cr3. 
for  all  Distinct  Element,  Finite  Element  and  Experimental  results:  (a)  in 
the  direction  of  <j\  and  (b)  in  the  direction  of  <r2  (a2  is  kept  constant). 
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Normalized  constrained  moduli,  versus  stress  ratio,  K  =  an  / 

<70,  for  the  case  of  biaxial  compression— extension  loading  with  constant 
mean  stress,  <r0  =  91  KPa. 


Figure  45.  The  RPI  Soil  Mechanics  Laboratory  MTS  Servohydraulic  Axial/ Torsional 
and  Computer  control  System  Configuration. 


Figure  46. 


Sketch  of  Hollow  Cylindrical  Specimen  Used  in  this  Study. 
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Figure  47. 


Grain  Size  Distribution  Curve  for  Mixture  of  Glass  Beads  Used  in  this 
Study. 
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Figure  49. 


Volumetric  Strain  versus  Octahedral  Shear  Strain  for  Monotonic  Radial 
Shear  Tests  (0  =  0°,  30°,  and  60°)  on  Assemblies  of  Different  Sizes  of 
Spheres. (crm  ~  140  KPa). 
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Figure  50.  q/p  Ratio  versus  Maximum  Shear  Strain  for  Monotonic  Radial  Shear  Tests 
Having  9  =  0°. 


octahedral  strain  (  *  ) 


Volumetric  Strain  versus  Octahedral  Shear  Strain  from  Monotonic 
Shear  Tests  ($  =  0°)  on  Mixtures  of  Different  Sizes  of  Glass 
Consolidated  at  Different  Confining  Stresses. 


octahedral  stress  (  KPa 


Sequences  of  Loading  for  the  Determination  of  the  "initial"  and 
"subsequent"  Yield  Surfaces. 


Figure  54. 


Flowchart  of  the  program  controlling  the  experiments 
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Figure  56.  Experimentally  Obtained  Initial  and  Subsequent  Yield  Surfaces  in  Glass 
Beads  (70Ct  =  0.05%). 
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Figure  58.  Initial  and  Subsequent  yield  surfaces  of  the  531— Sphere  Medium  in  the  rv h 
1/2  {o^-a  h)  Space  calculated  Using  Distinct  Element  Program 

CONBAL— 2.  Yield  Surfaces  Defined  as  the  Set  of  all  Points  with 

7oct  =  0.05%. 


Figure  60.  Initial  and  Subsequent  Yield  Surfaces  of  the  531-Sphere  Medium  in  the  rvh 
1/2  (av-0  h)  Space  calculated  Using  Distinct  Element  Program 

CONBAL— 2.  Yield  Surfaces  Defined  as  the  Locus  of  all  Points  with  70Ct  = 
0.05%. 


Figure  61.  Conical  Yield  Surfaces  and  Current  r— plane  for  Constitutive  Law  of 
Granular  Soil  of  Eqs.  1—7.  Position  of  Yield  Cones  Correspond  to  Loading 
from  0  to  A  followed  by  Stress  increment  AB.  Point  B  is  on  the  x=plane. 
The  Value  of  p  at  the  Apex  of  Any  Cone  (e.g.,  A,  B)  is  pi. 


(a) 


loading  path 


(b) 


Figure  62. 


Schematic  Representation  of  (a)  Translation  and  (b)  Distortion  of  a  Yield 
Surface  (Yen  1979). 


Figure  63. 


Predicted  and  Measured  Yield  Surfaces  in  Aluminum  (Yen  and  Eisenberg 
1987). 


Yield  Surface  in  principal  Stress  Space.  Note  the  Rotation  of  the  Surface 
Around  its  Apex  (Prevost  1985). 


Figure  64. 


